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We review the phenomenon of equilibrium fluctuations in the number of condensed atoms no in a 
trap containing N atoms total. We start with a history of the Bose-Einstein distribution, a similar 
grand canonical problem with an indefinite total number of particles, the Einstein-Uhlenbeck debate 
concerning the rounding of the mean number of condensed atoms no near a critical temperature 
T c , and a discussion of the relations between statistics of BEC fluctuations in the grand canonical, 
canonical, and microcanonical ensembles. 

First, we study BEC fluctuations in the ideal Bose gas in a trap and explain why the grand 
canonical description goes very wrong for all moments ((no — no)" 1 ), except of the mean value. 
We discuss different approaches capable of providing approximate analytical results and physical 
insight into this very complicated problem. In particular, we describe at length the master equation 
and canonical-ensemble quasiparticle approaches which give the most accurate and physically trans- 
parent picture of the BEC fluctuations. The master equation approach, that perfectly describes 
even the mesoscopic effects due to the finite number N of the atoms in the trap, is quite similar 
to the quantum theory of the laser. That is, we calculate a steady-state probability distribution 
of the number of condensed atoms p no {t = oo) from a dynamical master equation and thus get 
the moments of fluctuations. We present analytical formulas for the moments of the ground-state 
occupation fluctuations in the ideal Bose gas in the harmonic trap and arbitrary power-law traps. 

In the last part of the review, we include particle interaction via a generalized Bogoliubov for- 
malism and describe condensate fluctuations in the interacting Bose gas. In particular, we show 
that the canonical-ensemble quasiparticle approach works very well for the interacting gases and 
find analytical formulas for the characteristic function and all cumulants, i.e., all moments, of the 
condensate fluctuations. The surprising conclusion is that in most cases the ground-state occupa- 
tion fluctuations are anomalously large and are not Gaussian even in the thermodynamic limit. We 
also resolve the Giorgini, Pitaevskii and Stringari (GPS) vs. Idziaszek et al. debate on the vari- 
ance of the condensate fluctuations in the interacting gas in the thermodynamic limit in favor of 
GPS. Furthermore, we clarify a crossover between the ideal-gas and weakly- interacting-gas statistics 
which is governed by a pair-correlation, squeezing mechanism and show how, with an increase of 
the interaction strength, the fluctuations can now be understood as being essentially 1 /2 that of an 
ideal Bose gas. We also explain the crucial fact that the condensate fluctuations are governed by a 
singular contribution of the lowest energy quasiparticles. This is a sort of infrared anomaly which is 
universal for constrained systems below the critical temperature of a second order phase transition. 
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I. INTRODUCTION 



Professor Herbert Walther has taught us that good physics unifies and unites seemingly different fields. Nowhere is 
this more apparent than in the current studies of Bose-Einstein condensation (BEC) and coherent atom optics which 
draw from and contribute to the general subject of coherence effects in many-body physics and quantum optics. It 
is in this spirit that the present paper presents the recent application of techniques, ideas, and theorems which have 
been developed in understanding lasers and squeezed states to the condensation of N bosons. Highlights of these 
studies, and related points of BEC history, are described in the following paragraphs. 

1. Bose Jjl got the ball rolling by deriving the Planck distribution without using classical electrodynamics, 
as Planck [!| and Einstein Q had done. Instead, he took the extreme photon-as-a- particle point of view, and by 
regarding these particles as indistinguishable obtained, among other things, Planck's result, 



where is the mean number of photons with energy and wavevector k, [3 = (feaT) -1 , T is the blackbody 
temperature, and ks is Boltzmann's constant. 

However, his paper was rejected by the Philosophical Magazine and so he sent it to Einstein, who recognized its 
value. Einstein translated it into German and got it published in the Zeitschrift fur Physik Q]. He then applied Bose's 
method to atoms and predicted that the atoms would "condense" into the lowest energy level when the temperature 
was low enough 0, IE • 

Time has not dealt as kindly with Bose as did Einstein. As is often the case in the opening of a new field, things 
were presented and understood imperfectly at first. Indeed Bose did his "counting" of photon states in cells of phase 
space in an unorthodox fashion. So much so that the famous Max Delbruck wrote an interesting article @ in which 
he concluded that Bose made a mistake, and only got the Planck distribution by serendipity. We here discuss this 
opinion, and retrace the steps that led Bose to his result. Sure, he enjoyed a measure of luck, but his mathematics 
and his derivation were correct. 

2. Einstein's treatment of BEC of atoms in a large box showed a cusp in the number of atoms in the ground state, 
no, as a function of temperature, 



for T < T c , where N is the total number of atoms, and T c is the (critical) transition temperature. 

Uhlenbeck criticized this aspect of Einstein's work, claiming that the cusp at T = T c is unphysical. Einstein 
agreed with the Uhlenbeck criticism but argued that in the limit of large numbers of atoms (the thermodynamic limit) 
everything would be okay. Later, Uhlenbeck and his student Kahn showed |lfl ] that Einstein was right and put the 
matter to rest (for a while). 

Fast forward to the present era of mesoscopic BEC physics with only thousands (or even hundreds) of atoms in a 
condensate. What do we now do with this Uhlenbeck dilemma? As one of us (W. K.) showed some time ago all 
that is needed is a better treatment of the problem. Einstein took the chemical potential to be zero, which is correct 
for the ideal Bose gas in the thermodynamic limit. However, when the chemical potential is treated more carefully, 
the cusp goes away, as we discuss in detail, see, e.g., Figs.^, and0 

3. So far everything we have been talking about concerns the average number of particles in the condensate. Now 
we turn to the central focus of this review: fluctuations in the condensate particle number. As the reader will recall, 
Einstein used the fluctuation properties of waves and particles to great advantage. In particular he noted that in 
Planck's problem, there were particle-like fluctuations in photon number in addition to the wave-like contribution, 
i.e. 




(2) 




nji(wave) + (particle), 



(3) 
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and in this way he argued for a particle picture of light. 

In his studies on Bose-Einstein condensation, he reversed the logic arguing that the fluctuations in the ideal quantum 
gas also show both wave-like and particle-like attributes, just as in the case of photons. It is interesting that Einstein 
was led to the wave nature of matter by studying fluctuations. We note that he knew of and credited de Broglie at 
this point well before wave mechanics was developed. 

Another important contribution to the problem of BEC fluctuations came from Fritz London's observation that 
the specific heat is proportional to the variance of a Bose-Einstein condensate and showed a cusp, which he calculated 
as being around 3. IK. It is noteworthy that the so-called lambda point in liquid Helium, marking the transition from 
normal to superfluid, takes place at around 2.19K. 

However, Ziff, Uhlenbeck and Kac ^3 note several decades later that there is a problem with the usual treatment 
of fluctuations. They say: 

[When] the grand canonical properties for the ideal Bose gas are derived, it turns out that some of them 
differ from the corresponding canonical properties — even in the bulk limit! . . . The grand canonical 
ensemble . . . loses its validity for the ideal Bose gas in the condensed region. 

One of us (M. H.) has noted elsewhere |14j that: 

This grand canonical fluctuation catastrophe has been discussed by generations of physicists . . . 

Let us sharpen the preceding remarks. Large fluctuations are a feature of the thermal behavior of systems of bosons. 
If n is the mean number of non-interacting particles occupying a particular one particle state, then the mean square 
occupation fluctuation in the grand canonical picture is n(n + I). If, however, the system has a fixed total number of 
particles ./V confined in space by a trapping potential, then at low enough temperature T when a significant fraction 
of TV are in the ground state, such large fluctuations are impossible. No matter how large N, the grand canonical 
description cannot be even approximately true. This seems to be one of the most important examples that different 
statistical ensembles give agreement or disagreement in different regimes of temperatures. To avoid the catastrophe, 
the acclaimed statistical physicist D. ter Haar |l5j| proposed that the fluctuations in the condensate particle number 
in the low temperature regime (adapted to a harmonic trap) might go as 



(4) 



This had the correct zero limit as T — > 0, but is not right for higher temperatures where the leading term actually 
goes as [N (^-j^J ] 1 ^ 2 - The point is that fluctuations are subtle; even the ideal Bose gas is full of interesting physics 
in this regard. 

In this paper, we resolve the grand canonical fluctuation catastrophe in several ways. In particular, recent ap- 
plication of techniques developed in the quantum theory of the laser 0, 0] and in quantum optics [TsL ITflj allow 
us to formulate a consistent and physically appealing analytical picture of the condensate fluctuations in the ideal 
and interacting Bose gases. Our present understanding of the statistics of the BEC fluctuations goes far beyond the 
results that were formulated before the 90's BEC boom, as summarized by Ziff, Uhlenbeck, and Kac in their classical 
review [l^. Theoretical predictions for the BEC fluctuations, which are anomalously large and non-Gaussian even in 
the thermodynamical limit, are derived and explained on the basis of the simple analytical expressions |2Ctl2f| . The 
results are in excellent agreement with the exact numerical simulations. The existence of the infrared singularities in 
the moments of fluctuations and the universal fact that these singularities are responsible for the anomalously large 
fluctuations in BEC, are among the recent conceptual discoveries. The quantum theory of laser threshold behavior 
constitutes another important advance in the physics of bosonic systems. 

4. The laser made its appearance in the early 60's and provided us with a new source of light with a new kind of 
photon statistics. Before the laser, the statistics of radiation were either those of black-body photons associated with 
Planck's radiation, which for a single mode of frequency v takes the form 

Pn = e -nPHu( 1 _ e -PHu^ (5) 

or when one considers radiation from a coherent oscillating current such as a radio transmitter or a microwave klystron 
the photon distribution becomes Poissonian, 

77™ 

Pn = ^e"« (6) 
n! 

where n is the average photon number. 



5 




FIG. 1: (a) Mean value (no) and (b) variance Ano = \J (wf) — ( n o) 2 of the number of condensed atoms as a function of 
temperature for N = 200 atoms in a harmonic trap calculated via the solution of the condensate master equation (solid line). 
Large dots are the exact numerical results obtained in the canonical ensemble. Dashed line for (no) is a plot of N[l — (T/T c ) 3 ] 
which is valid in the thermodynamic limit. Dashed line for Ano is the grand canonical answer y/ no (no + 1) which gives 
catastrophically large fluctuations below T c . 



However, laser photon statistics, as derived from the quantum theory of the laser, goes from black-body statistics 
below threshold to Poissonian statistics far above threshold. In between, when we are in the threshold region (and 
even above threshold as in the case, for example, of the helium neon laser), we have a new distribution. We present 
a review of the laser photon statistical problem. 

It has been said that the Bose-Einstein condensate is to atoms what the laser is to photons; even the concept of 
an atom laser has emerged. In such a case, one naturally asks, "what is the statistical distribution of atoms in the 
condensate?" For example, let us first address the issue of an ideal gas of N atoms in contact with a reservoir at 
temperature T. The condensate occupation distribution in the harmonic trap under these conditions at low enough 
temperatures is given by the BEC master equation analysis as 

_ 1 [N(T/T c ) 3 ] N - no m 
Pno ~Z N (N-n )\ ■ [ '> 

The mean number and variance obtained from the condensate master equation are in excellent agreement with 
computer simulation (computer experiment) as shown in Fig.^ We will discuss this aspect of the fluctuation problem 
in some detail and indicate how the fluctuations change when we go to the case of the interacting Bose gas. 

5. The fascinating interface between superfluid He II and BEC in a dilute gas was mapped out by the experiments 
of Reppy and coworkers [5^: and finite-size effects were studied theoretically by M. Fisher and co-workers |23|. They 
carried out experiments in which He II was placed in a porous glass medium which serves to keep the atoms well 
separated. These experiments are characterized by a dilute gas BEC of N atoms at temperature T. 

Of course, it was the successful experimental demonstration of Bose-Einstein condensation in the ultracold atomic 
alkali-metal pil I2M |2(| , hydrogen j2?J and helium gases IH, that stimulated the renaissance in the theory of 
BEC. In less than a decade, many intriguing problems in the physics of BEC, that were not studied, or understood 
before the 90's [U EJ HUHl HI Hg , were formulated and resolved. 

6. Finally, we turn on the interaction between atoms in the BEC and find explicit expressions for the characteristic 
function and all cumulants of the probability distribution of the number of atoms in the (bare) ground state of a 
trap for the weakly interacting dilute Bose gas in equilibrium. The surprising result is that the BEC statistics is 
not Gaussian, i.e., the ratio of higher cumulants to an appropriate power of the variance does not vanish, even in 
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the thermodynamic limit. We calculate explicitly the effect of Bogoliubov coupling between excited atoms on the 
suppression of the BEC fluctuations in a box ( "homogeneous gas" ) at moderate temperatures and their enhancement 
at very low temperatures. We find that there is a strong pair-correlation effect in the occupation of the coupled 
atomic modes with the opposite wavevectors k and — k. This explains why the ground-state occupation fluctuations 
remain anomalously large to the same extent as in the non-interacting gas, except for a factor of 1/2 suppression. 
We find that, roughly speaking, this is so because the atoms are strongly coupled in correlated pairs such that the 
number of independent stochastic occupation variables ( "degrees of freedom" ) contributing to the fluctuations of 
the total number of excited atoms is only 1/2 the atom number N. This is a particular feature of the well studied 
quantum optics phenomenon of two-mode squeezing (see, e.g., [37j and 0,0]). The squeezing is due to the quantum 
correlations that build up in the bare excited modes via Bogoliubov coupling and is very similar to the noise squeezing 
in a nondegenerate parametric amplifier. 

Throughout the review, we will check main approximate analytical results (such as in Eqs. (|162[) . I|172|) . 1)2231) . (|263[) . 
(|271|l ) against the "exact" numerics based on the recursion relations i|79|) and l|80() which take into account exactly 
all mesoscopic effects near the critical temperature T c . Unfortunately, the recursion relations are known only for the 
ideal Bosc gas. In the present review we discuss the BEC fluctuations mainly in the canonical ensemble, which cures 
misleading predictions of the grand canonical ensemble and, at the same time, does not have any essential differences 
with the microcanonical ensemble for most physically interesting quantities and situations. Moreover, as we discuss 
below, the canonical partition function can be used for an accurate calculation of the microcanonical partition function 
via the saddle-point method. 



II. HISTORY OF THE BOSE-EINSTEIN DISTRIBUTION 



In late 1923, a certain Satyendranath Bose, reader in physics at the University of Dacca in East Bengal, submitted 
a paper on Planck's law of blackbody radiation to the Philosophical Magazine. Six months later he was informed 
that the paper had received a negative referee report, and consequently been rejected [2^. While present authors 
may find consolation in the thought that the rejection of a truly groundbreaking paper after an irresponsibly long 
refereeing process is not an invention of our times, few of their mistreated works will eventually meet with a recognition 
comparable to Bose's. Not without a palpable amount of self-confidence, Bose sent the rejected manuscript to Albert 
Einstein in Berlin, together with a handwritten cover letter dated June 4, 1924, beginning |39l |: 

Respected Sir: 

I have ventured to send you the accompanying article for your perusal and opinion. I am anxious to know 
what you think of it. You will see that I have tried to deduce the coefficient 8m/ 2 /c 3 in Planck's Law 
independent of the classical electrodynamics, only assuming that the ultimate elementary regions in the 
phase-space has the content h 3 . I do not know sufficient German to translate the paper. If you think the 
paper worth publication I shall be grateful if you arrange its publication in Zeitschrift fur Physik. Though 
a complete stranger to you, I do not hesitate in making such a request. Because we are all your pupils 
though profiting only from your teachings through your writings . . . 

In hindsight, it appears curious that Bose drew Einstein's attention only to his derivation of the prefactor in Planck's 
law. Wasn't he aware of the fact that his truly singular achievement, an insight not even spelled out explicitly in 
Einstein's translation of his paper as it was received by the Zeitschrift fur Physik on July 2, 1924 0,0, but contained 
implicitly in the mathematics, lay elsewhere? 



What Bose did 



In the opening paragraph of his paper Q, , Bose pounces on an issue which he considers unsatisfactory: When 
calculating the energy distribution of blackbody radiation according to 

8m/ 2 di/ , , 

Q v dv = - — E v , (8) 



that is, 



energy per volume of blackbody radiation with frequency between v and v + dv 
= number of modes contained in that frequency interval of the radiation field per volume 
x thermal energy E v of a radiation mode with frequency v , 
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the number of modes had previously been derived only with reference to classical physics. In his opinion, the logical 
foundation of such a recourse was not sufficiently secure, and he proposed an alternative derivation, based on the 
hypothesis of light quanta. 

Considering radiation inside some cavity with volume V, he observed that the squared momentum of such a light 
quantum is related to its frequency through 

where h denotes Planck's constant, and c is the velocity of light. Dividing the frequency axis into intervals of length 
dv s , such that the entire axis is covered when the label s varies from s = to s = oo, the phase space volume 
associated with frequencies between v and v + dv a therefore is 

dx dy dz dp x dpydp z — V in f — ^ 

h 3 v 2 

= VAir^dv 3 . (10) 

c 3 

It does not seem to have bothered Bose that the concept of phase space again brings classical mechanics into play. 
Relying on the assumption that a single quantum state occupies a cell of volume h 3 in phase space, a notion which, in 
the wake of the Bohr-Sommerfeld quantization rule, may have appeared natural to a physicist in the early 1920s, and 
accounting for the two states of polarization, the total number A s of quantum cells belonging to frequencies between 
v and v + dv s , corresponding to the number of radiation modes in that frequency interval, immediately follows: 

C 3 

That's all, as far as the first factor on the r.h.s. of Eq. © is concerned. This is what Bose announced in his letter to 
Einstein, but this is, most emphatically, not his main contribution towards the understanding of Planck's law. The 
few lines which granted him immortality follow when he turns to the second factor. Back-translated from Einstein's 
phrasing of his words 0, Q : 

Now it is a simple task to calculate the thermodynamic probability of a (macroscopically defined) state. 
Let iV s be the number of quanta belonging to the frequency interval dv s . How many ways are there to 
distribute them over the cells belonging to dv s l Let p s be the number of vacant cells, p\ the number of 
those containing one quantum, p| the number of cells which contain two quanta, and so on. The number 
of possible distributions then is 

A s| Sttv 2 

, where A s = V — —dv s , (12) 



pV-pV- 



and where 



A " • />;, • 1 • • 2 • ... 

is the number of quanta belonging to dv s . 

What is happening here? Bose is resorting to a fundamental principle of statistical mechanics, according to which 
the probability of observing a state with certain macroscopic properties — in short: a macrostate — is proportional 
to the number of its microscopic realizations — microstates — compatible with the macroscopically given restrictions. 
Let us, for example, consider a model phase space consisting of four cells only, and let there be four quanta. Let us 
then specify the macrostate by requiring that one cell remain empty, two cells contain one quantum each, and one cell 
be doubly occupied, i.e., po = 1, p\ = 2, pi — 1, and p r = for r > 3. How many microstates are compatible with 
this specification? We may place two quanta in one out of four cells, and then choose one out of the remaining three 
cells to be the empty one. After that there is no further choice left, since each of the two other cells now has to host 
one quantum. Hence, there are 4 x 3 = 12 possible configurations, or microstates: 12 = 4!/(l! 2! 1!). In general, when 
there are A s cells belonging to dv s , they can be arranged in A s \ ways. However, if a cell pattern is obtained from 
another one merely by a rearrangement of those p$ cells containing no quantum, the configuration remains unchanged. 
Obviously, there are Pq\ such "neutral" rearrangements which all correspond to the same configuration. The same 
argument then applies, for any r > 1, to those p s r cells containing r quanta: Each of the p s r \ possibilities of arranging 
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the cells with r quanta leads to the same configuration. Thus, each configuration is realized by p^.pl ! . . . equivalent 
arrangements of cells, and the number of different configurations, or microstates, is given by the total number of 
arrangements divided by the number of equivalent arrangements, that is, by Bose's expression Ijl2(l . 

There is one proposition tacitly made in this way of counting microstates which might even appear self-evident, 
but which actually constitutes the very core of Bose's breakthrough, and which deserves to be spelled out explicitly: 
When considering equivalent arrangements as representatives of merely one microstate, it is implied that the quanta 
are indistinguishable. It does not matter "which quantum occupies which cell"; all that matters are the occupation 
numbers p s r . Even more, the "which quantum" -question is rendered meaningless, since there is, as a matter of 
principle, no way of attaching some sort of label to individual quanta belonging to the same &v s , with the purpose 
of distinguishing them. This "indistinguishability in principle" does not occur in classical physics. Two classical 
particles may have the same mass, and identical other properties, but it is nevertheless taken for granted that one 
can tell one from the other. Not so, according to Bose, with light quanta. 

The rest of Bose's paper has become a standard exercise in statistical physics. Taking into account all frequency 
intervals d^ s , the total number of microstates corresponding to a pre-specified set {p s r } of cell occupation numbers is 

The logarithm of this functional yields the entropy associated with the considered set {pf,}. Since, according to the 
definition of p s r , 

A s = Y\* for each s , (14) 

r 

and assuming the statistically relevant p* to be large, Stirling's approximation Inn! ~ nlnn — n gives 

In W [{#}] = VJ A s In A* - ]T ]T p* hip* . (15) 

s s r 

The most probable macrostate now is the one with the maximum number of microstates, characterized by that set of 
occupation numbers which maximizes this expression lJT5)l. Stipulating that the radiation field be thermally isolated, 
so that its total energy 

E = H N ' f " / ' With N ' = X! ( 16 ) 

s r 

is fixed, the maximum is found by variation of the p s r , subject to this constraint (|16fl . In addition, the constraints 1141) 
have to be respected. Introducing Lagrangian multipliers X s for these "number-of-cells" constraints, and a further 
Lagrangian multiplier (3 for the energy constraint, the maximum is singled out by the condition 



5 (in W[{p S r }} " £ A' - /? £>' £ rtf) = 

\ s r ST/ 



(17) 



giving 



S P s r {bxp s r + l + X s ) + l3j2 h ^J2 r5p r = - ( 18 ) 



Since the Sp^ can now be taken as independent, the maximizing configuration {pf.} obeys 

lnp£ + 1 + X s + r{3hv s =0, (19) 

or 

= B s e- rf3hl,s , (20) 



with normalization constants B s to be determined from the constraints Q14JI: 
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The total number of quanta for the maximizing configuration then is 



A s 1 



E 



rc 



-rf3hv s 



A» 



(22) 



Still, the physical meaning of the Lagrangian multiplier (3 has to be established. This can be done with the help of 
the entropy functional, since inserting the maximizing configuration yields the thermodynamical equilibrium entropy: 



S = fc B lnW[{^}] 



(3E - AS ln f 1 - e 



-0hv' 



(23) 



where fee denotes Boltzmann's constant. From the identity dS/dE = 1/T one then finds (3 — l/(kBT), the inverse 
energy equivalent of the temperature T. Hence, from Eqs. I|22ll and Bose obtains the total energy of the radiation 
contained in the volume V in the form 



e = j2 NSh,yS 

S 

^ 8nhv s3 



ex P (B 



k-eT 



1) 



-Av s 



(24) 



which is equivalent to Planck's formula: With the indistinguishability of quanta, i.e., Bose's enumeration <|12ll of 
microstates as key input, the principles of statistical mechanics immediately yield the thermodynamic properties of 
radiation. 



B. What Einstein did 



Unlike that unfortunate referee of the Philosophical Magazine, Einstein immediately realized the power of Bose's 
approach. Estimating that it took the manuscript three weeks to travel from Dacca to Berlin, Einstein may have 
received it around June 25 Only one week later, on July 2, his translation of the manuscript was officiallyreceived 
by the Zeitschrift fur Physik. The author's name was lacking its initials — the byline of the published paper [lj simply 
reads: By Bose (Dacca University, India) — but otherwise Einstein was doing Bose fair justice: He even sent Bose 
a handwritten postcard stating that he regarded his paper as a most important contribution; that postcard seems 
to have impressed the German Consulate in Calcutta to the extent that Bose's visa was issued without requiring 
payment of the customary fee [38|. 

Within just a few days, Einstein then took a further step towards exploring the implications of the "indistinguish- 
ability in principle" of quantum mechanical entities. At the end of the printed, German version of Bose's paper 
there appears the parenthetical remark "Translated by A. Einstein", followed by an announcement: 

Note added by the translator: Bose's derivation of Planck's formula constitutes, in my opinion, an impor- 
tant step forward. The method used here also yields the quantum theory of the ideal gas, as I will explain 
in detail elsewhere. 

"Elsewhere" in this case meant the Proceedings of the Prussian Academy of Sciences. In the session of the Academy 
on July 10, Einstein delivered a paper entitled "Quantum theory of the monoatomic ideal gas" pj. In that paper, he 
considered nonrelativistic free particles of mass m, so that the energy-momentum relation simply reads 

E = ^> ^ 
and the phase-space volume for a particle with an energy not exceeding E is 

4-7T 

$ = U— {2mEf 2 . (26) 
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Again relying on the notion that a single quantum state occupies a cell of volume h 3 in phase space, the number of 
such cells belonging to the energy interval from E to E + AE is 

As = '^{2mf' 2 E 1 / 2 AE . (27) 

Thus, for particles with nonzero rest mass As is the analog of Bose's A s introduced in Eq. (|llfl . Einstein then specified 
the cell occupation numbers by requiring that, out of these As cells, p s r As cells contain r particles, so that p s r is the 
probability of finding r particles in any one of these cells, 



Now comes the decisive step. Without attempt of justification or even comment, Einstein adopts Bose's way <|12|) 
of counting the number of corresponding microstates. This is a far-reaching hypothesis, which implies that, unlike 
classical particles, atoms of the same species with energies in the same range AE are indistinguishable: Interchanging 
two such atoms does not yield a new microstate; as with photons, it does not matter "which atom occupies which 
cell". Consequently, the number of microstates associated with a pre-specified set of occupation probabilities {pf.} for 
the above As cells is 



giving, with the help of Stirling's formula, 



As! 

W ' = TT°° f SA M > 29 



In W = - As Pr ln Pr ■ (30) 



Einstein then casts this result into a more attractive form. Stipulating that the index s does no longer refer jointly to 
the cells within a certain energy interval, but rather labels individual cells, the above expression naturally generalizes 
to 

lnW{{p s r }} = -J2Y,Pr ln Pr> ( 31 ) 

s r 

where the cell index s now runs over all cells, so that p a r here is the probability of finding r particles in the s-th cell, 
ft is interesting to observe that this functional 131fl has precisely the same form as the Shannon entropy introduced 
in 1948 in an information-theoretical context [40j . 

Since ^2 r rp s T gives the expectation value of the number of particles occupying the cell labelled s, the total number 
of particles is 



while the total energy of the gas reads 



where E s is the energy of a particle in the s-th cell. Since, according to Eq. I|26|) . a cell's number s is related to the 
energy E s through 

S = V 

= ^(2min 3/2 , (34) 

one has 

E s = as 2 ' 3 (35) 

with 

-2/3 



h 2 / 4irV \ 
2^ \~~1T J 



(36) 
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Considering an isolated system, with given, fixed particle number TV and fixed energy E, the macrostate realized in 
nature is characterized by that set {p s r } which maximizes the entropy functional subject to the constraints H32|) 
and <|33[) . together with the constraints l|28|) expressing normalization of the cell occupation probabilities. Hence, 



\ s r s r s r / 



(37) 



so that, seen from the conceptual viewpoint, the only difference between Bose's variational problem (I17|) and Einstein's 
variational problem (|37() is the appearance of an additional Lagrangian multiplier a in the latter: In the case of 
radiation, the total number of light quanta adjusts itself in thermal equilibrium, instead of being fixed beforehand; 
in the case of a gas of particles with nonzero rest mass, the total number of particles is conserved, requiring the 
introduction of the entailing multiplier a. One then finds 

\np s r + 1 + A s + ar + f3rE s = 



or 



(38) 
(39) 
(40) 

Here we deviate from the; notation in Einstein's paper 5], in order to be compatible with modern conventions. The 
expectation value for the occupation number of the cell with energy E s then follows from an elementary calculation 
similar to Bose's reasoning Q22p: 



with normalization constants to be determined from the constraints J5SJ: 



B s 



-(a+(3E s ) 



1 



e a+f3E° _ i 



Therefore, the total number of particles and the total energy of the gas can be expressed as 

1 



N 



E 



E = E- 



e a+0E° 
E s 



e a+f3E* _ i 



(41) 

(42) 
(43) 



Inserting the maximizing set (|39|l into the functional (|31|l yields, after a brief calculation, the equilibrium entropy of 
the gas in the form 



S = k B lnW[{p s r }} 



aN - 



PE-^bxfl-e-^ 



ht3E s 



(44) 



In order to identify the Lagrangian multiplier /3, Einstein considered an infinitesimal heating of the system, assuming 
its volume and, hence, the cell energies E s to remain fixed. This gives 



dE 



TdS 



Nda + pdE + Edp~Y,%^ 



k B Tf3dE 



requiring 



k B T 



(45) 



(46) 



As in Bose's case, the Lagrangian multiplier (5 accounting for the energy constraint is the inverse energy equivalent of 
the temperature T. The other multiplier a, guaranteeing particle number conservation, then is determined from the 
identity 

In the following two sections of his paper Q, Einstein shows how the thermodynamics of the classical ideal gas 
is recovered if one neglects unity against e a+ P E , and derives the virial expansion of the equation of state for the 
quantum gas obeying Eqs. I|42|) and (|43|l . 
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C. Was Bose— Einstein statistics arrived at by serendipity? 



The title of this subsection is a literal quote from the title of a paper by M. Delbriick [8j, who contents that Bose 
made an elementary mistake in statistics in that he should have bothered "which quantum occupies which cell" , which 
would have been the natural approach, and that Einstein first copied that mistake without paying much attention 
to it. Indeed, such a suspicion does not seem to be entirely unfounded. In his letter to Einstein, Bose announces 
only his comparatively straightforward derivation of the number (| 1 II) of radiation modes falling into the frequency 
range from v to v + d^ s , apparently being unaware that his revolutionary deed was the implicit exploitation of the 
"indistinguishability in principle" of quanta — a concept so far unheard of. In Einstein's translation of his paper pj 
this notion of indistinguishability does not appear in words, although it is what underlies the breakthrough. Even 
more, it does not appear in the first paper Q on the ideal Bose gas — until the very last paragraph, where Einstein 
ponders over 

...a paradox which I have been unable to resolve. There is no difficulty in treating also the case of a 
mixture of two different gases by the method explained here. In this case, each molecular species has its 
own "cells" . From this follows the additivity of the entropies of the mixture's components. Therefore, with 
respect to molecular energy, pressure, and statistical distribution each component behaves as if it were the 
only one present. A mixture containing n\ and molecules, with the molecules of the first kind being 
distinguishable (in particular with respect to the molecular masses mi, m^) only by an arbitrarily small 
amount from that of the second, therefore yields, at a given temperature, a pressure and a distribution of 
states which differs from that of a uniform gas with ri\ + ri2 molecules with practically the same molecular 
mass, occupying the same volume. However, this appears to be as good as impossible. 

Interestingly, Einstein here considers "distinguishability to some variable degree" , which can be continuously reduced 
to indistinguishability. But this notion is flawed: Either the molecules have some feature which allows us to tell one 
species from the other, in which case the different species can be distinguished, or they have none at all, in which case 
they are indistinguishable in principle. Thus, at this point, about two weeks after the receipt of Bose's manuscript 
and one week after sending its translation to the Zeitschrift fur Physik, even Einstein may not yet have fully grasped 
the implications of Bose's way (11211 of counting microstates. 

But there was more to come. In December 1924, Einstein submitted a second manuscript on the quantum theory of 
the ideal Bose gas 0,0, formally written as a continuation of the first one. He began that second paper by pointing 
out a curiosity implied by his equation of state of the ideal quantum gas: Given a certain number of particles N and 
a temperature T, and considering a compression of the volume V, there is a certain volume below which a segregation 
sets in. With decreasing volume, an increasing number of particles has to occupy the first quantum cell, i.e., the state 
without kinetic energy, while the rest is distributed over the other cells according to Eq. I|41|l . with e° = 1. Thus, 
Bose-Einstein condensation was unveiled! But this discovery merely appears as a small addendum to the previous 
paper Q, for Einstein then takes up a different, more fundamental scent. He mentioned that Ehrenfest and other 
colleagues of his had criticized that in Bose's and his own theory the quanta or particles had not been treated as 
statistically independent entities, a fact which had not been properly emphasized. Einstein agrees, and then he sets 
out to put things straight. He abandons his previous "single-cell" approach and again considers the collection of 
quantum cells with energies between E v and E v + AE V , the number of which is 

Zv = ™(2mf' 2 El/^E u . (47) 

Then he juxtaposes in detail Bose's way of counting microstates to what is done in classical statistics. Assuming that 
there are n v quantum particles falling into Ai?„, Bose's approach l|12|) implies that there are 

W V = ^ + (48) 
nj [z v - 1)! 

possibilities of distributing the particles over the cells. This expression can easily be visualized: Drawing the n v 
particles as a sequence of n v "dots" in a row, they can be organized into a microstate with specific occupation 
numbers for z v cells — again assuming that it does not matter which particle occupies which cell — by inserting 
z v — 1 separating "lines" between them. Thus, there are n v + z v — 1 positions carrying a symbol, n v of which are dots. 
The total number of microstates then equals the total number of possibilities to select the n v positions carrying a 
"dot" out of these n v + z v — 1 positions, which is just the binomial coefficient stated in Eq. Q48)). To give an example: 
Assuming that there are n„ = 4 particles and z v = 4 cells, Eq. (|48|l states that there are altogether 
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microstates. On the other hand, there are several sets of occupation numbers which allow one to distribute the 
particles over the cells: 



occupation numbers 


number of microstates 


P4 = 1, po = 3 


1!3I — 4 


P3 = 1, Pi = 1, Po = 2 


4! - 12 

1!1!2! - iz 


P2 = 2, Po = 2 


-JL = 6 

2! 2! U 


P2 = 1, Pi = 2, J5 = 1 


4! - 12 

1!2M! - iz 


pi = 4 


il - 1 

4! - 1 



The right column of this table gives, for each set, the number of microstates according to Bose's formula (|12H ; 
obviously, these numbers add up to the total number 35 anticipated above. Thus, the binomial coefficient 
conveniently accounts for all possible microstates, without the need to specify the occupation numbers, according to 
the combinatorial identity 



E 



z\ 



Po! 



Pat! 



N + Z -1 
N 



(49) 



where the sum is restricted to those sets {po,Pi, ■ ■ ■ ,Pn} which comply with the two conditions ^2 r p r = Z and 
y] rp r — N, as in the example above. In Appendix^ we provide a proof of this identity. With this background, let 
us return to Einstein's reasoning: When taking into account all energy intervals AE U , the total number of microstates 
is given by the product W = Y[ v W v , providing the entropy functional 



lnW[{n„}] = y] \(n v + z v ) ln(n v + z v ) — n v In n v — z v In z v 



E 







+ z v In ( — - 


-)] 




n u J 







The maximizing set {n v } now has to obey the two constraints 

J2 n - = N 

n v E v = E , 



(50) 

(51) 
(52) 



but there is no more need for the multipliers A s appearing in the previous Eqs. i|17|) and Q37[l. since the constraints l|14|) 
or (|28|l are automatically respected when starting from the convenient expression l|48|) . Hence, one has 



5 ^ln W[{n u }\ -aJ2 n u~P^2 n » E A = 



(53) 



leading immediately to 



E 



In 1 



— a — {3E V 



5n v = , 



e a+/3E v _ I 



(54) 



(55) 



in agreement with the previous result (|41|l . 

But what, Einstein asks, would have resulted had one not adopted Bose's prescription <|12t and thus counted 
equivalent arrangements with equal population numbers only once, but rather had treated the particles as classical, 
statistically independent entities? Then there obviously are 



W v = (z v ) 



(56) 
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possibilities of distributing the n v particles belonging to AE V over the z v cells: Each particle simply is placed in one 
of the z v cells, regardless of the others. Now, when considering all intervals AE V , with distinguishable particles it 
does matter how those n v particles going into the respective AE V are selected from all N particles; for this selection, 
there are N\/Y\ v n v l possibilities. Thus, taking classical statistics seriously, there are 



W = AH J| 



(57) 



possible microstatcs, yielding 



lnW[{rc„}] = NlnN - N + [n v \nz v 



n v In n v + n u 





+ n v 







(58) 



This is a truly vexating expression, since it gives a thermodynamical entropy which is not proportional to the total 
number of particles, i.e., no extensive quantity, because of the first term on the r.h.s. Hence, already in the days before 
Bose and Einstein one had got used to ignoring the leading factor Nl in Eq. (|57|l . with the half-hearted concession 
that microstates which result from each other by a mere permutation of the N particles should not be counted as 
different. Of course, this is an intrinsic inconsistency of the classical theory: Instead of accepting that, shouldn't one 
abandon Eq. I|57|l straight away and accept the more systematic quantum approach, despite the apparently strange 
consequence of losing the particles' independence? And Einstein gives a further, strong argument in favor of the 
quantum theory: At zero temperature, all particles occupy the lowest cell, giving n± = N and n v = for v > 1. With 
z\ = 1, the quantum way of counting based on Eq. I|48|> gives just one single microstate, which means zero entropy in 
agreement with Nernst's theorem, whereas the classical expression (|57f) yields an incorrect entropy even if one ignores 
the disturbing AM Finally, the variational calculation based on the classical functional 1)58(1 proceeds via 



E 



In I — )- a-/3E v 

71) 



5n, y 







(59) 



furnishing the Boltzmann-like distribution 



= z„e- a -P E » 



(60) 



for the maximizing set {n,}. In short, the quantum ideal gas of nonzero mass particles with its distribution 1551) 
deviates from the classical ideal gas in the same manner as does Planck's law of radiation from Wien's law. This 
observation convinced Einstein, even in the lack of any clear experimental evidence, that Bose's way of counting 
microstates had to be taken seriously, since, as he remarks in the introduction to his second paper "if it is justified 
to consider radiation as a quantum gas, the analogy between the quantum gas and the particle gas has to be a complete 
one" . This belief also enabled him to accept the sacrifice of the statistical independence of quantum particles implied 
by the formula 148(1 . which, by the end of 1924, he had clearly realized: 

The formula therefore indirectly expresses a certain hypothesis about a mutual influence of the molecules 
on each other which is of an entirely mysterious kind ... 

But what might be the physics behind that mysterious influence which non-interacting particles appear to exert 
on each other? In a further section of his paper 0, Einstein's reasoning takes an amazing direction: He considers 
the density fluctuations of the ideal quantum gas, and from this deduces the necessity to invoke wave mechanics! 
Whereas he had previously employed what is nowadays known as the microcanonical ensemble, formally embodied 
through the constraints that the total number of particles and the total energy be fixed, he now resorts to a grand 
canonical framework and considers a gas within some finite volume V which communicates with a gas of the same 
species contained in an infinitely large volume. He then stipulates that both volumes be separated from each other 
by some kind of membrane which can be penetrated only by particles with an energy in a certain infinitesimal range 
AE V , and quantifies the ensuing fluctuation An v of the number of particles in V, not admitting energy exchange 
between particles in different energy intervals. Writing n v = n u + An„, the entropy of the gas within V is expanded 
in the form 



Sg as {An u ) — Sg 



dSg as 
dAn v 



An, 



1 d 2 S gas 

2 d{An v f 



(An,) 2 



(61) 
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whereas the entropy of the reservoir changes with the transferred particles according to 

Sb(An„) =S - T^An„ . (62) 
c/An, 

In view of the assumed infinite size of the reservoir, the quadratic term is negligible here. Since the equilibrium state 
is characterized by the requirement that the total entropy S = Sg^s + Sq be maximum, one has 

i,S , (63) 



SAn, 

so that 

1 d 2 'S, 



S(An v ) = S + £ " ~ ga ° (An,) 2 . (64) 



2 d(An,) 2 

Hence, the probability distribution for finding a certain fluctuation An, is Gaussian 

P(An„) = const ■ e S(A "" )/fcB 



2fc B c>(An,) 2 

from which one reads off the mean square of the fluctuations, 

\2\ _ fc B 



cons t ■ exp ( -L f / Sl " (An,) 2 ) , (65) 



((An,) 2 ) = fl ° ^ . (66) 



Since, according to the previous Eq. I|54[l. one has 



1 dSg as 



/cb <9An, \ n 

one deduces 

1 9 2 S 



'3(An„) 



= ln(l + ^l , (67) 



^gas 



fcs 3(An,) 2 n 2 + z,n, 

resulting in 



(68) 



<(A7^) 2 ) = n, + ^ (69) 

or 

((An,/n,) 2 ) = 1- + - . (70) 
n, z, 

With z v = 1, this gives the familiar grand canonical expression for the fluctuation of the occupation number of a 
single quantum state. With a stroke of genius, Einstein now interprets this formula: Whereas the first term on the 
r.h.s. of Eq. (|70|l would also be present if the particles were statistically independent, the second term reminds him 
of interference fluctuations of a radiation field 6]: 

One can interpret it even in the case of a gas in a corresponding manner, by associating with the gas a 
radiation process in a suitable manner, and calculating its interference fluctuations. I will explain this in 
more detail, since I believe that this is more than a formal analogy. 

He then refers to de Broglie's idea of associating a wavelike process with single material particles and argues that, if 
one associates a scalar wave field with a gas of quantum particles, the term 1 / z v in Eq. Q70[l describes the corresponding 
mean square fluctuation of the wave field. What an imagination — on the basis of the fluctuation formula l|7U|l Einstein 
anticipates many-body matter waves, long before wave mechanics was officially enthroned! Indeed, it was this paper 
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of his which led to a decisive turn of events: From these speculations on the relevance of matter waves Schrodinger 
learned about de Broglie's thesis, acquired a copy of it, and then formulated his wave mechanics. 

Having recapitulated this history, let us once again turn to the title of this subsection: Was Bose-Einstein statistics 
arrived at by serendipity? Delbriick's opinion that it arose out of an elementary mistake in statistics that Bose made 
almost certainly is too harsh. On the other hand, the important relation l|49(l . see also the Appendix lAl does not 
seem to have figured in Bose's thinking. When writing down the crucial expression l|12(l . Bose definitely must have 
been aware that he was counting the number of microstates by determining the number of different distributions of 
quanta over the available quantum cells, regardless of "which quantum occupies which cell" . He may well have been 
fully aware that his way of counting implied the indistinguishability of quanta occupying the same energy range, but 
he did not reflect on this curious issue. On the other hand, he didn't have to, since his way of counting directly led 
to one of the most important formulas in physics, and therefore simply had to be correct. 

Many years later, Bose recalled |41j: 

I had no idea that what I had done was really novel .... I was not a statistician to the extent of really 
knowing that I was doing something which was really different from what Boltzmann would have done, 
from Boltzmann statistics. Instead of thinking of the light-quantum just as a particle, I talked about these 
states. 

By counting the number of ways to fill a number of photonic states (cells) Bose obtained Eq. (|13|) which is exactly 
the same form as Boltzmann's Eq. <|57[) for z v = 1, but with new meanings attached to the new symbols: n s replaced 
by p s r and N replaced by A s . Bose's formula leads to an entirely different new statistics - the quantum statistics for 
indistinguishable particles, in contrast to Boltzmann's distinguishable particles. It took a while for this to sink in 
even with Einstein, but that is the nature of research. 

Contrary to Bose, Einstein had no experimental motivation when adapting Bose's work to particles with nonzero 
rest mass. He seems to have been guided by a deep-rooted belief in the essential simplicity of physics, so that he was 
quite ready to accept a complete analogy between the gas of light quanta and the ideal gas of quantum particles, 
although he may not yet have seen the revolutionary implications of this concept when he submitted his first paper p| 
on this matter. But the arrival at a deep truth on the basis of a well-reflected conviction can hardly be called 
serendipity. His second paper [|J is, by all means, a singular intellectual achievement, combining daring intuition 
with almost prophetical insight. And who would blame Einstein for trying to apply, in another section of that second 
paper, his quantum theory of the ideal Bose gas to the electron gas in metals? 

In view of the outstanding importance which Einstein's fluctuation formula (|70|l has had for the becoming of wave 
mechanics, it appears remarkable that a puzzling question has long remained unanswered: What happens if one faces, 
unlike Einstein in his derivation of this relation 170|) . a closed system of Bose particles which does not communicate 
with some sort of particle reservoir? When the temperature T approaches zero, all N particles are forced into the 
system's ground state, so that the mean square ((Ano) 2 ) of the fluctuation of the ground-state occupation number 
has to vanish for T — > — but with zq = 1, the grand canonical Eq. (|69|l gives ((Ano) 2 ) - * N(N + 1), clearly 
indicating that with respect to these fluctuations the different statistical ensembles are no longer equivalent. What, 
then, would be the correct expression for the fluctuation of the ground-state occupation number within the canonical 
ensemble, which excludes any exchange of particles with the environment, but still allows for the exchange of energy? 
To what extent does the microcanonical ensemble, which applies when even the energy is kept constant, differ from the 
canonical one? Various aspects of this riddle have appeared in the literature over the years (l3l IT5I l42| . mainly inspired 
by academic curiosity, before it resurfaced in 1996 |43l l44l l45l l46j . this time triggered by the experimental realization 
of mesoscopic Bose-Einstein condensates in isolated microtraps. Since then, much insight into this surprisingly rich 
problem has been gained, and some of the answers to the above questions have been given by now. In the following 
sections of this article, these new developments will be reviewed in detail. 



D. Comparison between Bose's and Einstein's counting of the number of microstates W 

In Bose's original counting 1|12|) . he considered the numbers pf, of cells occupied with r photons, so that the total 
number of cells is given by A s = Y^Pr- While Einstein still had adopted this way of counting in his first paper [I| 

r 

on the ideal Bose gas, as witnessed by his Eq. I|29[) . he used, without comment, more economical means in his second 
paper J6j, relying on the binomial coefficient l|48l) . Figure[2]shows an example in which A^ s — 2 particles are distributed 
over A s = 3 cells, and visualizes that Bose's formula for counting the number of possible arrangements (or the number 
of microstates which give the same macrostate) gives the same result as Einstein's formula only after one has summed 
over all possible configurations, i.e., over those setspg,pf , ... which simultaneously obey the two conditions A s = J^r-Pr 
and A^ s = ^2 r rp*. In this example, there are only two such sets; a general proof is given in Appendix lAl Thus, it 
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degeneracy (number of states) , g=3 
number of particles , n=2 



p =2, p^O, p 2 =l p =l,P!=2,p 2 =0 

3 ! 3! 

W Bose j = - = 3 w Bose 2 = = 



2!0!1! 1!2!0! 

\YBose= -\YBose i+ W Bose 2 = 3+3= 6 

(g+n-1)! 4! 

W Einstein = = =6 = W Bose 

n!(g-l)! 2!2! 

FIG. 2: A simple example showing Bose's and Einstein's (textbook) methods of counting the number of possible ways to put 
n = 2 particles in a level having A — g — 3 states. 

is actually possible to state the number of microstates without specifying the individual arrangements, by summing 
over all of them: 

y ... y -r-P — , A° = y P ° and jv* = y>s] = (71) 

^ ^ p^-pV- ■ ■ ■ Pir, ^ ^ N S UA S -1)\ v ' 



III. GRAND CANONICAL VERSUS CANONICAL STATISTICS OF BEC FLUCTUATIONS 

The problem of BEC fluctuations in a Bose gas is well known [l3l . ll5| . However, as with many other problems which 
are well-posed physically and mathematically, it is highly nontrivial and deep, especially for the interacting Bose gas. 

A. Relations between statistics of BEC fluctuations in the grand canonical, canonical, and microcanonical 

ensembles 

To set the stage, we briefly review here some basic notions and facts from the statistical physics of BEC involving 
relations between statistics of BEC fluctuations in the grand canonical, canonical, and microcanonical ensembles. 

Specific experimental conditions determine which statistics should be applied to describe a particular system. In 
view of the present experimental status the canonical and microcanonical descriptions of the BEC are of primary 
importance. Recent BEC experiments on harmonically trapped atoms of dilute gases deal with a finite and well 
defined number of particles. This number, even if it is not known exactly, certainly does not fluctuate once the 
cooling process is over. Magnetic or optical confinement suggests that the system is also thermally isolated and, 
hence, the theory of the trapped condensate based on a microcanonical ensemble is needed. In this microcanonical 
ensemble the total particle number and the total energy are both exactly conserved, i.e., the corresponding operators 
are constrained to be the c-number constants, N = const and H = const. On the other hand, in experiments with 
two (or many) component BECs, Bose-Fermi mixtures, and additional gas components, e.g., for sympathetic cooling, 
there is an energy exchange between the components. As a result, each of the components can be described by the 
canonical ensemble that applies to systems with conserved particle number while exchanging energy with a heat bath 
of a given temperature. Such a description is also appropriate for dilute 4 He in a porous medium [22j. In the canonical 
ensemble only the total number of particles is constrained to be an exact, non- fluctuating constant, 

N = h + J2h k , (72) 

but the energy H has non-zero fluctuations, ((H — H) 2 )) ^ 0, and only its average value is a constant, (H) = H — 
E = const, determined by a fixed temperature T of the system. 
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However, the microcanonical and canonical descriptions of a many-body system are difficult because of the operator 
constraints imposed on the total energy and particle number. As a result, the standard textbook formulation of the 
BEC problem assumes either a grand canonical ensemble (describing a system which is allowed to exchange both energy 
and particles with a reservoir at a given temperature T and chemical potential /i, which fix only the average total 
number of particles, (N) — const, and the average total energy (H ) = const) or some restricted ensemble that selects 
states so as to ensure a condensate wave function with an almost fixed phase and amplitude [3fl ItH l4^ | . These 
standard formulations focus on and provide effective tools for the study of the thermodynamic and hydrodynamic 
properties of the many-body Bose system at the expense of an artificial modification of the condensate statistics 
and dynamics of BEC formation. While the textbook grand canonical prediction of the condensate mean occupation 
agrees, in some sense, with the Bose-Einstein condensation of trapped atomic gases, this is not even approximately 
true as concerns the grand canonical counting statistics 

P? C M = -^(-%-) n \ (73) 



1 + flu V 1 

which gives the probability to find n v particles in a given single-particle state v, where the mean occupation is n v . 
Below the Bosc-Einstein condensation temperature, where the ground state mean-occupation number is macroscopic, 
fio ~ N, the distribution pg (no) becomes extremely broad with the squared variance ((no — no) 2 ) ~ fig oc N 2 even at 
T — > [l3ll49l |. This prediction is surely at odds with the isolated Bose gas, where for sufficiently low temperature all 
particles are expected to occupy the ground state with no fluctuations left. It was argued by Ziff et al. that this 
unphysical behavior of the variance is just a mathematical artefact of the standard grand canonical ensemble, which 
becomes unphysical below the condensation point. Thus, the grand canonical ensemble is irrelevant to experiment if 
not revised properly. Another extremity, namely, a complete fixation of the amplitude and phase of the condensate 
wave function, is unable to address the condensate formation and the fluctuation problems at all. 

It was first realized by Fierz |42| that the canonical ensemble with an exactly fixed total number of particles removes 
the pathologies of the grand canonical ensemble. He exploited the fact that the description of BEC in a homogeneous 
ideal Bose gas is exactly equivalent to the spherical model of statistical physics, and that the condensate serves as 
a particle reservoir for the non-condensed phase. Recently BEC fluctuations were studied by a number of authors 
in different statistical ensembles, both in the ideal and interacting Bose gases. The microcanonical treatment of the 
ground-state fluctuations in a one-dimensional harmonic trap is closely related to the combinatorics of partitioning 
an integer, opening up an interesting link to number theory [43l l50l l5T | . 

It is worthwhile to compare the counting statistics (|73[) with the predictions of other statistical ensembles [5| . For 
high temperatures, T > T c , all three ensembles predict the same counting statistics. However, this is not the case 
for low temperatures, T < T c . Here the broad distribution of the grand canonical statistics differs most dramatically 
from po(^o) i n the canonical and microcanonical statistics, which show a narrow single-peaked distribution around 
the condensate mean occupation number. In particular, the master equation approach |52l l53| (Section IV) yields a 
finite negative binomial distribution for the probability distribution of the ground-state occupation in the ideal Bose 
gas in the canonical ensemble (see Fig. 1121 and Eqs. ^ 154(1 . (|162ll . and 11631) ). The width of the peak decreases with 
decreasing temperature. In fact, it is this sharply-peaked statistical distribution which one would naively expect for 
a Bose condensate. 

Each statistical ensemble is described by a different partition function. The microcanonical partition function 
Q(E, N) is equal to the number of TV-particle microstates corresponding to a given total energy E. Interestingly, at 
low temperatures canonical and microcanonical fluctuations have been found to agree in the large- N (thermodynamic) 
limit for one-dimensional harmonic trapping potentials, but to differ in the case of three-dimensional isotropic harmonic 
traps [54j]. More precisely, for the d-dimensional power-law traps characterized by an exponent tr, as considered 
later in Section V.C, microcanonical and canonical fluctuations agree in the lar ge- N limit when d/a < 2, but the 
microcanonical fluctuations remain smaller than the canonical ones when d/a > 2 [55(. Thus, in d = 3 dimensions the 
homogeneous Bose gas falls into the first category, but the harmonically trapped one into the second. The difference 
between the fluctuations in these two ensembles is expressed by a formula which is similar in spirit to the one expressing 
the familiar difference between the heat capacities at constant pressure and at constant volume |54tl55|. 

The direct numerical computation of the microcanonical partition function becomes very time consuming or not 
possible for N > 10 5 . For N 3> 1, and for large numbers of occupied excited energy levels, one can invoke, e.g., 
the approximate technique based on the saddle-point method, which is widely used in statistical physics |3l|. When 
employing this method, one starts from the known grand canonical partition function and utilizes the saddle-point 
approximation for extracting its required canonical and microcanonical counterparts, which then yield all desired 
quantities by taking suitable derivatives. Recently, still another statistical ensemble, the so called Maxwell's demon 
ensemble, has been introduced |54| . Here, the system is divided into the condensate and the particles occupying excited 
states, so that only particle transfer (without energy exchange) between these two subsystems is allowed, an idea that 
had also been exploited by Fierz |42| and by Politzer |44j. This ensemble has been used to obtain an approximate 
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analytical expression for the ground-state BEC fluctuations both in the canonical and in the microcanonical ensemble. 
The Maxwell's demon approximation can be understood on the basis of the canonical-ensemble quasiparticle approach 
[2(1 l2lj | , which is discussed in Section V, and which is readily generalizable to the case of the interacting Bose gas 
(see Section VI). It also directly proves that the higher statistical moments for a homogeneous Bose gas depend on 
the particular boundary conditions imposed, even in the thermodynamic limit [2Cll2lLl56| . 
The canonical partition function Zn{T) is defined as 

oo 

Z N {T) = Y / e- E/kBT n{E,N), (74) 

E 

where the sum runs over all allowed energies, and k,B is the Boltzmann constant. Eq. (|74() can be used to calculate also 
the microcanonical partition function Q(E, N) by means of the inversion of this definition. Likewise, the canonical 
partition function Zn (T) can be obtained by the inversion of the definition of the grand canonical partition function 
S( M ,T): 

OO 

S(ji,T)= J2^ N/kBT Z N (T). (75) 

JV=0 

Inserting Eq. (|74|l into Eq. I|75|) . we obtain the following relation between S(/i, T) and fl(E,N): 

oo oo 

EOh T)=J2 ^ N/ksT J2 e~ E,kBT n(E, N). (76) 

JV=0 E 

As an example, let us consider an isotropic 3-dimensional harmonic trap with cigcnfrcqucncy lu. In this case, a 
relatively compact expression for the grand canonical partition function, 

oo , . s (E/Hu+l)(E/Hul+2)/2 

allows us to find fi(_E, N) from Eq. H7f>|) by an application of the saddle-point approximation to the contour integral 

"(S. *0 = fojjjg £ dz £ dx Z ^/Li • x = e-*"»*, z = (78) 

where the contours of integration 7 2 and ^ x include z — and x — 0, respectively. It is convenient to rewrite the 
function under the integral in Eq. (|78Jl as exp[<£>(z,:c)], where (p(z,x) — lnS(z,x) — (N + l)lnz — [E/fyjj + l)hix. 
Taking the contours through the extrema (saddle points) zq and xo of ip{z,x), and employing the usual Gaus- 
sian approximation, we get for N — > oo and E/fku — * oo the following asymptotic formula: Q(E,N) = 

[2ttD(zq, xo)] _1//2 S(zo, %o) I ' Zn +1 x E ^ hu ' +1 , where D(zq, xq) is the determinant of the second derivatives of the function 
ip(z,x), evaluated at the saddle points |45|. For N — > oo and Ejfvjj — > oo the function exp[<p(z, x)] is sharply peaked 
at zq and xq, which ensures good accuracy of the Gaussian approximation. However^ the standard result becomes 
incorrect for E < Ne±, where e\ = fiwi is the energy of the first excited state (see |i^,|53); in this case, a more re- 
finded version of the saddle-point method is required [5§L l59l l60| . We discuss this improved version of the saddle-point 
method in Appendix IF1 

An accurate knowledge of the canonical partition function is helpful for the calculation of the microcanonical 
condensate fluctuations by the saddle point method, as has been demonstrated [57| by a numerical comparison 
with exact microcanonical simulations. In principle, the knowledge of the canonical partition function allows us to 
calculate thermodynamic and statistical equilibrium properties of the system in the standard way (see, e.g., [T^.l6l|). 
An important fact is that the usual thermodynamic quantities (average energy, work, pressure, heat capacities, etc.) 
and the average number of condensed atoms do not have any infrared-singular contributions and do not depend on 
a choice of the statistical ensemble in the thermodynamic (bulk) limit. However, the variance and higher moments 
of the BEC fluctuations do have the infrared anomalies and do depend on the statistical ensemble, so that their 
calculation is much more involved and subtle. 



B. Exact recursion relation for the statistics of the number of condensed atoms in an ideal Bose gas 



It is worth noting that there is one useful reference result in the theory of BEC fluctuations, namely, an exact 
recursion relation for the statistics of the number of condensed atoms in an ideal Bose gas. Although it does not give 
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any simple analytical answer or physical insight into the problem, it can be used for "exact" numerical simulations 
for traps containing a finite number of atoms. This is very useful as a reference to be compared against different 
approximate analytical formulas. This exact recursion relation for the ideal Bose gas had been known for a long 
time and re-derived independently by several authors [46l l62l l63l l64| . In the canonical ensemble, the probability 
to find uq particles occupying the single-particle ground state is given by 

, \ Z N-n Q {T) - Z N _ no _i(T) 

Po{n ) = T~FF\ ; Z -! = ■ ( 79 ) 

Zn{-1 ) 



The recurrence relation for the ideal Bose gas then states |61l |6 

N 

X 



-. N oo 

ZNi^^-^ZiiT/^Z^T), Z l {T)=Y j e- e " ,kBT , Z Q (T) = 1, (80) 



k=l u=Q 



which enables one to numerically compute the entire counting statistics (|79H . Here v stands for a set of quantum 
numbers which label a given single-particle state, e v is the associated energy, and the ground-state energy is taken as 
£o = by convention. For an isotropic, three-dimensional harmonic trap one has 



Z 1 {T) = Y j \{n + 2){n + l)e-^ 



ihuj/ksT _ f 

where \[n + 2)(n + 1) is the degeneracy of the level e n — nhcu. 

In the microcanonical ensemble the ground-state occupation probability is given by a similar formula 

Po (no)- , O(£?,-l) = 0, 



where the microcanonical partition function obeys the recurrence relation 

N oo 

X 



N oo 

n(E,N) = — ^2^Tn(E- ks v ,N- k), n(0,N) = l, n(E>0,0) = 0. (82) 



fe=l v=0 

For finite E the sum over v is finite because of H(E < 0, N) = 0. 



C. Grand canonical approach 

Here we discuss the grand canonical ensemble, and show that it loses its validity for the ideal Bose gas in the 
condensed region. Nevertheless, reasonable approximate results can be obtained if the canonical-ensemble constraint 
is properly incorporated in the grand canonical approach, especially if we are not too close to T c . In principle, the 
statistical properties of BECs can be probed with light (23. In particular, the variance of the number of scattered 
photons may distinguish between the Poisson and microcanonical statistics. 

1. Mean number of condensed particles in an ideal Bose gas 



The Bose-Einstein distribution can be easily derived from the density matrix approach. Consider a collection of 

k 

system is described by 



particles with the Hamiltonian H = aj,afc(£fc — where fi is the chemical potential. The equilibrium state of the 



p=±exp(-0fi) (83) 

where Z =Tr{exp(— pH)} = Yiki^ ~ e~^ £fc ~''- ) ) _1 . Then the mean number of particles with energy is 

(n k ) = Tr{paia fe } = (l-e-^-^)^(n fe |a t fe a fe e-' 3a ^^-' i )|n fe ) 

s d(l - p-^Ofc-A'h-i 1 

= (l- e -^ ek -^)— - — = - (84) 

1 ' d(-P(e k -fx)) exp[/3( £fe -M)]-l' 1 ' 
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In the grand canonical ensemble the average condensate particle number no is determined from the equation for 
the total number of particles in the trap, 

oo oo 

-»)-! ■ < 85 > 

where is the energy spectrum of the trap. In particular, for the three dimensional (3D) isotropic harmonic trap we 
have £k = hfl(k x + k y + k z ). For simplicity, we set e = 0. 

For 3D and ID traps with non-interacting atoms, Eq. I|85|) was studied by Ketterle and van Druten and 
by Grossmann and Holthaus [43l l66l| . They calculated the fraction of ground-state atoms versus temperature for 
various N and found that BEC also exists in ID traps, where the condensation phenomenon looks very similar to 
the 3D case. Later Herzog and Olshanii used the known analytical formula for the canonical partition function 
of bosons trapped in a ID harmonic potential [6^, |6^| and showed that the discrepancy between the grand canonical 
and the canonical predictions for the ID condensate fraction becomes less than a few per cent for N > 10 4 . The 
deviation decreases according to a 1/lniV scaling law for fixed T/T c . In 3D the discrepancy is even less than in the 
ID system |64j . The ground state occupation number and other thermodynamic properties were studied by Chase, 
Mckjian and Zamick |6J] in the grand canonical, canonical and microcanonical ensembles by applying combinatorial 
techniques developed earlier in statistical nuclear fragmentation models. In such models there are also constraints, 
namely the conservation of proton and neutron number. The specific heat and the occupation of the ground state 
were found substantially in agreement in all three ensembles. This confirms the essential validity of the use of the 
different ensembles even for small groups of particles as long as the usual thermodynamic quantities, which do not 
have any infrared singular contributions, are calculated. 

Let us demonstrate how to calculate the mean number of condensed particles for a particular example of a 3D 
isotropic harmonic trap. Following Eq. I|85[l . we can relate the chemical potential fi to the mean number of condensed 
particles hq as 1 + I /no — exp(— /3/i). Thus, we have 

iV = V(n k )=V 7 — . (86) 

ti ^ (l + l/n )exp/3 £k -l > 

The standard approach is to consider the case N 3> 1 and separate off the ground state so that Eq. (|86|l approximately 
yields 

For an isotropic harmonic trap with frequency f2, 

v 1 = i v (»+ 2 )(" +i ) ~ i r (z+2)(.x+i) f88) 

f-> exp(/3e fc ) - 1 2^ exp(/3nM7) -1^ 2j l exp(a^fifi) - 1 [ ' 

rC^^O Tit — 1 

In the limit feT> Ml we obtain 



E 



i i r°° x 2 fk R T^ 3 



o exp(/3£fc) — 1 2 Jo exp(a;/?fif2) — 1 \ HQ 
Furthermore, when T — T c we take no = 0. Then Eqs. ((HZJ and ^ yield 

N 



-£7T C(3). (89) 



k B T c - nn ^— j (90) 

and the temperature dependence of the mean condensate occupation with a cusp at T = T c in the thermodynamic 
limit 




(91) 

Figure |3| compares the numerical solution of Eq. (|86|l (solid line) for N = 200 with the numerical calculation of 
fio(T) from the exact recursion relations in Eqs. 179f) and l|80l) in the canonical ensemble (large dots). Small dots 
show the plot of the solution l|91|l . which is valid only for a large number of particles, N. Obviously, more accurate 
solution of the equation for the mean number of condensed particles (I86|l in a trap with a finite number of particles 
does not show the cusp. In Appendix 151 we derive an analytical solution of Eq. (|86|l for fio(T) valid for fio(T) ^> 1. 
One can see that for the average particle number both ensembles yield very close answers. However, this is not the 
case for the BEC fluctuations. 
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FIG. 3: Mean number of condensed particles as a function of temperature for N = 200. The solid line is the plot of the 
numerical solution of Eq. I|86jl . The exact numerical result for the canonical ensemble (Eqs. I|79^ and JHOJ) is plotted as large 
dots. The dashed line is obtained using Maxwell's demon ensemble approach, which yields no = N — J2T>o V[ ex P(/3£fc) — !]■ 

2. Condensate fluctuations in an ideal Bose gas 

Condensate fluctuations are characterized by the central moments ((no — (no)) s ) = Yl r (^) ( n o)( n o) s_r - The first 
few of them are 

{(n -n ) 2 ) = (n 2 )-(n ) 2 , (92) 



((no - no) 3 ) = (ng) - 3<ng)<n<>) + 2(n ) 3 , (93) 



((no - n ) 4 ) - (n 4 ) - 4(ng) (no) + 6<ng) (n ) 2 - 3(n ) 4 . (94) 

At arbitrary temperatures, BEC fluctuations in the canonical ensemble can be described via a stochastic variable 
n = N — X^fc^o Uk * na * depends on and is complementary to the sum of the independently fluctuating numbers rife, 
k 7^ 0, of the excited atoms. In essence, the canonical ensemble constraint in Eq. I|72l) eliminates one degree of 
freedom (namely, the ground-state one) from the set of all independent degrees of freedom of the original grand 
canonical ensemble, so that only transitions between the ground and excited states remain independently fluctuating 
quantities. They describe the canonical-ensemble quasiparticles, or excitations, via the creation and annihilation 
operators j3 + and /3 (see Sections V and VI below). 

At temperatures higher than T c the condensate fraction is small and one can approximately treat the condensate 
as being in contact with a reservoir of non-condensate particles. The condensate exchanges particles with the big 
reservoir. Hence, the description of particle fluctuations in the grand canonical picture, assuming that the number of 
atoms in the ground state fluctuates independently from the numbers of excited atoms, is accurate in this temperature 
range. 

At temperatures close to or below T c the situation becomes completely different. One can say that at low tempera- 
tures, T <C T c , the picture is opposite to the picture of the Bose gas above the BEC phase transition at T > T c . The 
canonical-ensemble quasiparticle approach, suggested in Refs. ppl l2l| and valid both for the ideal and interacting 
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Bose gases, states that at low temperatures the non-condensate particles can be treated as being in contact with the 
big reservoir of the n condensate particles. This idea had previously been spelled out by Fierz J42| and been used by 
Politzer , and was then employed for the construction of the "Maxwell's demon ensemble" |54j . named so since a 
permanent selection of the excited (moving) atoms from the ground state (static) atoms is a problem resembling the 
famous Maxwell's Demon problem in statistical physics. 

Note that although this novel statistical concept can be studied approximately with the help of the Bose-Einstein 
expression (|84|l for the mean number of excited (only excited, k ^ 0!) states with some new chemical potential u, it 
describes the canonical-ensemble fluctuations and is essentially different from the standard grand canonical descrip- 
tion of fluctuations in the Bose gas. Moreover, it is approximately valid only if we are not too close to the critical 
temperature T c , since otherwise the "particle reservoir" is emptied. Also, the fluctuations obtained from the outlined 
"grand" canonical approach (see Eqs. I|95l) - (|102|l and l|C10|) - (|Cll|l '). an approach complementary to the grand canon- 
ical one, provide an accurate description in this temperature range in the thermodynamic limit, but do not take into 
account all mesoscopic effects, especially near T c (for more details, see Sections IV and V). Thus, although the mean 
number of condensed atoms no can be found from the grand canonical expression used in Eq. (|85(l . we still need to 
invoke the conservation of the total particle number N = no + n in order to find the higher moments of condensate 
fluctuations, (uq). Namely, we can use the following relation between the central moments of the roth order of the 
number of condensed atoms, and that of the non-condensed ones: ((no — no)" 1 ) — (— l) m ((n — n) m ). As a result, at 
low enough temperatures one can approximately write the central moments in the well-known canonical form via the 
cumulants in the ideal Bose gas (see Eqs. (|220|) . I|221l) . (|219|l . (|223|l and Section V below for more details): 



((n - n )' 



K2 = K2 



Kl 



k>0 



(n 2 k + n k ) 



(95) 



((n - n ) 3 ) = -k 3 = -(k 3 - 3k 2 - «i) 



fc>0 



(2n| + ml + ft fc ), 



(96) 



((n — no) 4 ) = ^4 + 3^2 = «4 + 6k 3 + 7k 2 + ki + 3(k 2 + «i) 



fc>0 



(6n^ + 12n^ + 7n^ + n fc ) + 3 



.fc>0 



(n 2 k + n k ) 



(97) 



These same equations can also be derived by means of the straightforward calculation explained in Appendix IU1 

Combining the hallmarks of the grand canonical approach, namely, the value of the chemical potential u — 
— ln(l + 1/n-o) and the mean non-condensate occupation n k — {exp[f3(Ek — /i)] — l} -1 , with the Eq. I|95|l describing 
the fluctuations in the canonical ensemble, we find the BEC variance 



An 2 = ((n - n ) 2 ) - £ 



fc>0 







1 




n j 





ex p(/3£ fe )(l + x) 



- f 



(98) 



In the case UbT 3> htt, this Eq. 198JI can be evaluated analytically, as is shown in Appendix El The variance up to 
second order in 1/no from Eq. 1D4|I is 



An 



N (T 



C(3) \T t 



G 
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n 



(1 + lnn ) + 









nl 




i)] 



-,1/2 



The leading term in this expression yields Politzer's result 44], 

An ? 



(99) 



(100) 



plotted as a dashed line in Fig.BJ where C(2) = y~l. 644 9 and C(3) ~ 1.202 1 (compare this with D. ter Haar's pj| 

result Ano ~ N i which is missing the square root). The same formula was obtained later by Navez et al. using 

the Maxwell's demon ensemble jS^]. For the microcanonical ensemble the Maxwell's demon approach yields [54[ 



(101) 
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FIG. 4: Variance ((no — (no)) 2 ) 1//2 of the condensate particle number for N = 200. The solid line is the "grand" canonical result 
obtained from Eq. 1981 and the numerical solution of (no) from Eq. 1861 . Exact numerical data for the canonical ensemble 
(Eqs. J23 and E3) 

are shown as dots. The asymptotic Politzer approximation j44|, given by Eq. 1 1 IKH . is shown by the dashed 
line, while small dots result from Eq. (ID4L The dash-dotted line is obtained from the master equation approach (see Eq. 1 17211 
below) . 



higher order terms have been derived in Ref. |60|. The microcanonical fluctuations are smaller than the canonical 
ones due to the additional energy constraint. For 2D and ID harmonic traps the Maxwell's demon approach leads to 

An ~ ^Ay^hr^ and An ~ , (102) 

respectively. 

Figure0]shows the BEC variance An (T) as a function of temperature for a 3D trap with the total particle number 
N = 200. The "grand" canonical curve refers to Eq. (|98(l and shows good agreement for T < T c with the numerical 
result for Ano(T) (large dots), obtained within the exact recursion relations ()79|) and l|80|l for the canonical ensemble. 
The plot of Politzer's asymptotic formula (jl(K)fl (dashed line) does not give good agreement, since the particle number 
considered here is fairly low. We also plot y/fho(na + 1), which is the expression for the condensate number fluctuations 
Ano in the grand canonical ensemble; it works well above T c . Figure [5] shows the third central moment ((no — no) 3 } 
as a function of temperature for the total particle number N = 200, plotted using Eqs. I|9rj|l and (|86[) . Dots are the 
exact numerical result obtained within the canonical ensemble. We also plot the standard grand canonical formula 
2ti,q + 3n,Q + no, which again works well only above T c . 

At high temperatures the main point, of course, is the validity of the standard grand canonical approach where 
the average occupation of the ground state alone gives a correct description of the condensate fluctuations, since 
the excited particles constitute a valid "particle reservoir". Condensate fluctuations obtained from the "grand" 
canonical approach for the canonical ensemble quasiparticles and from the standard grand canonical ensemble provide 
an accurate description of the canonical-ensemble fluctuations at temperatures not too close to the narrow crossover 
region between low (Eq. I|98|0 and high (a/ n,o{n,o + 1)) temperature regimes; i.e., in the region not too near T c . In 
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FIG. 5: The third central moment ((no — no) 3 ) for N = 200, plotted using Eqs. (I96H and 186H . Exact numerical data for the 
canonical ensemble (Eqs. 1791 and I8UI 1 are shown as dots. The dashed line is the result of master equation approach (see 
Eq. itTTH below). 



this crossover range both approximations fail, since the condensate and the non-condensate fractions have comparable 
numbers of particles, and there is no any valid particle reservoir. Note, however, that a better description, that 
includes mesoscopic effects and works in the whole temperature range, can be obtained using the condensate master 
equation, as shown in Section IV; see, e.g., Fig. 1131 Another (semi-analytic) technique, the saddle-point method, is 
discussed in Appendix IFl 



IV. DYNAMICAL MASTER EQUATION APPROACH AND LASER PHASE-TRANSITION ANALOGY 

One approach to the canonical statistics of ideal Bose gases, presented in |§2 an d developed further in J^, consists 
in setting up a master equation for the condensate and finding its equilibrium solution. This approach has the merit 
of being technically lucid and physically illuminating. Furthermore, it reveals important parallels to the quantum 
theory of the laser. In deriving that master equation, the approximation of detailed balance in the excited states is 
used, in addition to the assumption that given an arbitrary number no of atoms in the condensate, the remaining 
N — tiq excited atoms are in an equilibrium state at the prescribed temperature T. 

In Section IV. B we summarize the master equation approach against the results provided by independent techniques. 
In Section IV. A we motivate our approach by sketching the quantum theory of the laser with special emphasis on the 
points relevant to BEC. 



A. Quantum theory of the laser 

In order to set the stage for the derivation of the BEC master equation, let us remind ourselves of the structure of 
the master equations for a few basic systems that have some connection with N particles undergoing Bosc-Einstein 
condensation while exchanging energy with a thermal reservoir. 
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1. Single mode thermal field 

The dissipative dynamics of a small system (s) coupled to a large reservoir (r) is described by the reduced density 
matrix equation up to second order in the interaction Hamiltonian V sr 

^Ps(t) = -^2 Tr rj Q \yar(t),\y„tf),p.(t) (103) 

where /5*' 1 is the density matrix of the thermal reservoir, and we take the Markovian approximation. 

We consider the system as a single radiation mode cavity held (/) of frequency v coupled to a reservoir (r) of 
two-level thermal atoms, and show how the radiation held thcrmalizcs. The multiatom Hamiltonian in the interaction 
picture is 

V fr = h^gji&jcJe 1 ^^ + adj.) (104) 

3 

where dj = \bj){a,j \ is the atomic (spin) operator of the j-th particle corresponding to the downward transition, a* is 
the creation operator for the single mode held and gj is the coupling constant. The reduced density matrix equation 
for the held can be obtained from Eqs. 1)103(1 and 1)104(1 . By using the density matrix for the thermal ensemble of 
atoms 

pf = nd^X^"^ + \bj)(bj\e-^)/Z J (105) 

3 

where Zj = e ~ l3Ea ] + e _/3 - EbJ , we obtain 
d 1 

— pf = — - ^2 Kj{P aj (acft pf - 2a t pfa + pfaa*) + P bj (aJapf - 2apfd? + pfdJa)} (106) 

3 

where P Xj = e~^ Ex ^ /Zj with x — a.b and the dissipative constant is \nj —Re{gj f Q e^ -1 ^* - *'^'}. Note that the 
same structure of the master equation is obtained for a phonon bath modelled as a collection of harmonic oscillators, 
as shown in Appendix lEl 

Taking the diagonal matrix elements p n . n (t) — (n\pf\n) of Eq. 1(106(1 . we have 

d 

TT:Pn,n(<) = -nP a {{n + l)/9„, n (t) - up„_i, n _i(t)} - KP b {np n . n (t) - (n + n+1 (t)} , (107) 

at 

where k is Kj times a density of states factor and p n ,n- 
The steady state equation gives 

Pn,n=Pn = e- n0hu p O . (108) 

oo 

From ^2 Pn — 1 we obtain the thermal photon number distribution 

n=0 

p n = e - n ^{\ - e-^) (109) 
which is clearly an exponentially decaying photon number distribution. 



2. Coherent state 

Consider the interaction of a single mode held with a classical current J described by 

V coh (t) = [ J(r,t)-A(r,i)d 3 r 
Jv 

= h{j(t)a + f(t)tf) (110) 

where the complex time dependent coefficient is j(t) = ^ J v 3(r,t)-xe l ^ ]l ' r ^^ d 3 r and A a is the amplitude vector 
potential A(r,i) of the single mode held, assumed to be polarized along x axis. An example of such interaction is a 
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FIG. 6: Photon number distributions for a) thermal photons plotted from Eq. I|l()8fl (dashed line), b) coherent state (Poissonian) 
(thin solid line), and c) He-Ne laser plotted using Eq. I|121fl (thick solid line). Insert shows an atom making a radiation transition. 



klystron. Clearly, the unitary time evolution of V co h is in the form of a displacement operator exp(a* (t)d — a(t)a)) 
associated with a coherent state when dissipation is neglected. 

Thus, the density matrix equation for a klystron including coupling with a thermal bath is 



dp fit) 
dt 



1 



1 



= ^coh(i), Pf{t)] - -Tr r I [V fr (t), [V fr (t'), p f (t) <g> pf^dt' 



(111) 



where Vf r is given by Eq. (|E3|) . The second term of Eq. (|lll|l describes the damping of the single mode field given 
by Eq. (|E4|) . By taking the matrix elements of Eq. (|lllf> . after a bit of analysis, we have 



dp n ,n'(t) 
dt 



-«0'( < )v / n+ lp n +l,n' +j*(t)Vnp r , 



j(t)Vn'p 



lp. 



'n,n' + l 



1 



1 



--^■[(n + n')p ntn i - 1\J (n + l)(n' + l)p„+i,„'+i] - -V[(n + 1 + n' + l)p n ,n' - 2v / nn'Pn-i,n'-i] 



(112) 



Clearly, the first line of Eq. I|112|) shows that the change in the photon number is effected by the off-diagonal field 
density matrix or the coherence between two states of different photon number. On the other hand, the damping 
mechanism only causes a change in the photon number through the diagonal matrix element or the population of the 
number state. This is depicted in Fig. 

It can be shown that the solution of Eq. i|112[l for T> = is the matrix element of a coherent state \(3) = 

^ /2 E^W,i.e. 



Pn,n>(t) = (n\{\(3)(P\}\n'} = MWVr' e-W)? 

where f3(t) = a(t) — \C f Q a(t')dt' . This can be verified if we differentiate Eq. i|113|) 

dp n , n > f da(t) 1 r /— da*(t) 1 

—jj— = {— t -Ca(t)}Vnp n -i tn > + vn'{— t -Ca (i)}p n , n /_i 



(113) 



and using 



, daft) 

where ^ t 



~da(t) 
' dt 



-Ca(t)}Vn' + lp 



Vnpn 



n,n' + 1 



Pn-1, 



.da*(t) 
■ dt 



-Ca* (t)}Vn + Tp n+1 



Tl Pn,n' — Pn,n' — 1; 



Vn + lp n+ i^ n ' + i — p n ,n' + l, \frJ~+ lp„+l, rl ' + l — Pn+l,n> 

-ij*(t) is found by comparing Eq. (|114|l with Eq. 1)112(1 . 



(114) 



(115) 
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FIG. 7: Diagonal (star) and off-diagonal (circle) density matrix elements that govern temporal dynamics in a) klystron and b) 
thermal field. 




FIG. 8: Typical setup of a laser showing an ensemble of atoms driving a single mode field. A competition between lasing and 
dissipation through cavity walls leads to a detailed balance. 



Laser master equation 



The photon number equation i|l(J7|) for thermal field is linear in photon number, n, and it describes only the thermal 
damping and pumping due to the presence of a thermal reservoir. Now, we introduce a laser pumping scheme to drive 
the single mode field and show how the atom-field nonlinearity comes into the laser master equation. 

We consider a simple three level system where the cavity field couples level a and level b of lasing atoms in a 
molecular beam injected into a cavity in the excited state at rate r (see Fig. |SJ. The atoms undergo decay from level 
b to level c. The pumping mechanism from level c up to level a can thus produce gain in the single mode field. As 
shown in [13, we find 



dp n ,n(t) 
dt 



= -{A(n + 1) - B(n + l) 2 }p n ,n + {An - Bn 2 }p n -x,„-i, 



(116) 
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where A = ^f- is the linear gain coefficient and B — ^-A is the self-saturation coefficient. Here g is the atom- field 
coupling constant and 7 is the b — > c decay rate. We take the damping of the field to be 

d pnsi(t) \ _ _ Cnf)nn + + n+1 _ (117) 

°^ / loss 

Thus, the overall master equation for the laser is 
dp n ,n(t) 



dt 



= -{A(n + 1) - B(n + 1) }p n ,n + - Bn }p n -i,n-i ~ Cnp n , n + C(n + l)p n +i, n +i, (118) 



which is valid for small Bj A <C 1. 

We emphasize that the nonlinear process associated with B is a key physical process in the laser (but not in a 
thermal field) because the laser field is so large. 

We proceed with detailed balance equation between level n — 1 and n 

~{A(n + l)-B{n+l) 2 }p n +C(n+l) Pn+1 = 0, (119) 
{An - Bn 2 }p n -! - Cnp n = 0. (120) 

By iteration of p^ = ^ c Bk Pk-i, we have 

Pn=Pof[^^ (121) 

fc=l 

oo n oo 

where p = 1/(1 + II ) follows from P« = 1- Ec L- <|121[) is plotted in Fig. There we clearly see that 

n— 1 fc— 1 n— 

the photon statistics of, e.g. a He-Ne, laser is not Poissonian p n = (n) n e~' n ' /n\, as would be expected for a coherent 
state. 



B. Laser phase-transition analogy 

Bose- Einstein condensation of atoms in a trap has intriguing similarities with the threshold behavior of a laser 
which also can be viewed as a kind of a phase transition |7Ct l71| . Spontaneous formation of a long range coherent 
order parameter, i.e., macroscopic wave function, in the course of BEC second order phase transition is similar to 
spontaneous generation of a macroscopic coherent field in the laser cavity in the course of lasing. In both systems 
stimulated processes are responsible for the appearance of the macroscopic order parameter. The main difference is 
that for the Bose gas in a trap there is also interaction between the atoms which is responsible for some processes, 
including stimulated effects in BEC. Whereas for the laser there are two subsystems, namely the laser field and the 
active atomic medium. The crucial point for lasing is the interaction between the field and the atomic medium which 
is relatively small and can be treated perturbatively. Thus, the effects of different interactions in the laser system are 
easy to trace and relate to the observable characteristics of the system. This is not the case in BEC and it is more 
difficult to separate different effects. 

As is outlined in the previous subsection, in the quantum theory of laser, the dynamics of laser light is conveniently 
described by a master equation obtained by treating the atomic (gain) media and cavity dissipation (loss) as reservoirs 
which when "traced over" yield the coarse grained equation of motion for the reduced density matrix for laser radiation. 
In this way we arrive at the equation of motion for the probability of having n photons in the cavity given by Eq. (|118fl . 
From Eq. I|121|l we have the important result that partially coherent laser light has a sharp photon distribution (with 
width several times Poissonian for a typical He-Ne laser) due to the presence of the saturation nonlinearity, B, in the 
laser master equation. Thus, we see that the saturation nonlinearity in the radiation-matter interaction is essential 
for laser coherence. 

One naturally asks: is the corresponding nonlinearity in BEC due to atom-atom scattering? Or is there a nonlin- 
earity present even in an ideal Bose gas? The master equation presented in this Section proves that the latter is the 
case. 

More generally we pose the question: Is there a similar nonequilibrium approach for BEC in a dilute atomic gas 
that helps us in understanding the underlying physical mechanisms for the condensation, the critical behavior, and 
the associated nonlinearities? The answer to this question is "yes" ■ 
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FIG. 9: Simple harmonic oscillators as a thermal reservoir for the ideal Bose gas in a trap. 



C. Derivation of the condensate master equation 



We consider the usual model of a dilute gas of Bose atoms wherein interatomic scattering is neglected. This ideal 
Bose gas is confined inside a trap, so that the number of atoms, N, is fixed but the total energy, E, of the gas is not 
fixed. Instead, the Bose atoms exchange energy with a reservoir which has a fixed temperature T. As we shall see, 
this canonical ensemble approach is a useful tool in studying the current laser cooled dilute gas BEC experiments 
0IH El IHEi]. It is also directly relevant to the He-in-vycor BEC experiments ^ . 

This "ideal gas + reservoir" model allows us to demonstrate most clearly the master equation approach to the 
analysis of dynamics and statistics of BEC, and in particular, the advantages and mathematical tools of the method. 
Its extension for the case of an interacting gas which includes usual many-body effects due to interatomic scattering 
will be discussed elsewhere. 

Thus, we are following the so-called canonical-ensemble approach. It describes, of course, an intermediate situation 
as compared with the microcanonical ensemble and the grand canonical ensemble. In the microcanonical ensemble, 
the gas is completely isolated, E — const, N = const, so that there is no exchange of energy or atoms with a reservoir. 
In the grand canonical ensemble, only the average energy per atom, i.e., the temperature T and the average number 
of atoms (N) are fixed. In such a case there is an exchange of both energy and atoms with the reservoir. 

The "ideal gas + thermal reservoir" model provides the simplest description of many qualitative and, in some 
cases, quantitative characteristics of the experimental BEC. In particular, it explains many features of the condensate 
dynamics and fluctuations and allows us to obtain, for the first time, the atomic statistics of the BEC as discussed in 
the introduction and in the following. 



1, The "ideal gas + thermal reservoir" model 



For many problems a concrete realization of the reservoir system is not very important if its energy spectrum is 
dense and flat enough. For example, one expects (and we find) that the equilibrium (steady state) properties of the 
BEC are largely independent of the details of the reservoir. For the sake of simplicity, we assume that the reservoir is 
an ensemble of simple harmonic oscillators whose spectrum is dense and smooth, see Fig. ED. The interaction between 
the gas and the reservoir is described by the interaction picture Hamiltonian 

v = j2J29^ b M e ~ i(UJi ~ Uk+ul)t + H - c -> ( 122 ) 

j k>l 

where bj is the creation operator for the reservoir j oscillator ("phonon"), and a k and [k ^ 0) are the creation and 
annihilation operators for the Bose gas atoms in the fcth level. Here Hvu is the energy of the feth level of the trap, and 
gj t ki is the coupling strength. 



2. Bose gas master equation 



The motion of the total "gas + reservoir" system is governed by the equation for the total density matrix in the 
interaction representation, ptotal(t) = —i[V (t) , ptotalif)]/ h- Integrating the above equation for ptotah inserting it back 
into the commutator in Eq. (|123[) . and tracing over the reservoir, we obtain the exact equation of motion for the 
density matrix of the Bose-gas subsystem 

= -^i ^'TWm [V{t'),p total {t')]\, (123) 
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where Tr res stands for the trace over the reservoir degrees of freedom. 

We assume that the reservoir is large and remains unchanged during the interaction with the dynamical subsystem 
(Bose gas). As discussed in the density operator for the total system "gas+reservoir" can then be factored, i.e., 
Ptotai{t') ~ p(t') ® Pres, where p res is the equilibrium density matrix of the reservoir. If the spectrum is smooth, we 
are justified in making the Markov approximation, viz. p(f) —* p(t)- We then obtain the following equation for the 
reduced density operator of the Bose-gas subsystem, 

P = -^'^{ri k i + l)[a\aia\a k p~2a\a k pa\ai+ pa\aia\a k ] 
k>l 

-'^'^■n kl [a k a\a l a\p - 2aia\pa k a\ + pa k a\aia\]. (124) 

k>l 

In deriving Eq. (|124(l . we replaced the summation over reservoir modes by an integration (with the density of reservoir 
modes D(oj k i)) and neglected the frequency dependence of the coefficient n = 2irDg 2 /ft 2 . Here 

r\u = r](u kl ) = Tr res tf(uj k i)b(uj k i) = [exp(ftw fci /T) - 1] _1 (125) 

is the average occupation number of the heat bath oscillators at frequency u> k i = v k — v\ . Equation l|124|) is then the 
equation of motion for an N atom Bose gas driven by a heat bath at temperature T. 



3. Condensate master equation 

What we are most interested in is the probability distribution p no — J2{n k } Pn ,{n k }„ f° r the number of 
condensed atoms no, i.e., the number of atoms in the ground level of the trap. Let us introduce Pn ,{n k } nQ — 
(no, {n k }„ \p\no, {n k } no ) as a diagonal element of the density matrix in the canonical ensemble where + Sfc>o nk ~ 
N and \no,{n k } no ) is an arbitrary state of N atoms with occupation numbers of the trap's energy levels, n k , subject 
to the condition that there are no atoms in the ground state of the trap. 

In order to get an equation of motion for the condensate probability distribution p no , we need to perform the 
summation over all possible occupations {n k } no of the excited levels in the trap. The resulting equation of motion 
for p no , from Eq. (|124fl . is 

dpn 



dt 



{n k }„ k>l>0 

-ni{n k + 1K 0:{ . ..,„,_!,. ..,„ fe+ i : ... } „J 

+Vkl[ni(n k + l)Pn ,{n k } no 

-(ni + l)nkPn ,{...,m+i,...,n k -i,...} no ]} 

-« Y^k' + l ){ n + 1 ) n k>Pn .,{n k } no 

{n k } no fc'>0 

-(Vk> + l>o(n fc / + iK-i.in.+s^,}^., 

+m'M n k' + i)p„,,{ nt }„ 

-rjk'{n + l)n k 'P no+ i,{n h -5 h:k ,} na+1 ]> ( 126 ) 

where rj k > — r\[y k ^ is the mean number of thermal phonons of mode k' and the sum Y) k i runs over all excited levels. 

To simplify Eq. 11261) we assume that the atoms in the excited levels with a given number of condensed atoms no 
are in an equilibrium state at the temperature T, i.e., 

_ n exp(-A^ fc>0 ^n fc ) 

2^{„'J„ ex Pl T 2^k>0 v k n k) 

where X)fc>o nk = ^ ~ n 0i an d we assume that the sum ~^2 k>0 runs over all energy states of the trap, including 
degenerate states whose occupations n k are treated as different stochastic variables. Equation 1|127|) implies that the 
sum X)fc>z>o m Eq- l|126(l is equal to zero, since as depicted in Fig. 
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FIG. 10: Detailed balance and the corresponding probability flow diagram. We call K nQ cooling (kooling) rate since the laser 
loss rate is denote by "C". 



(Vkl + 1 )Pn ,{ni,}„ D - VklPn ,{...,m+l r 
{Vkl + l)Pn ,{...,n,-l,...,n fc +l,...} no = VklPn ,{n k } no ■ 



.,n k -l,— )„„ > 



(128) 



Equation l|128|) is precisely the detailed balance condition. The average number of atoms in an excited level, subject 
to the condition that there are uq atoms in the ground state, from Eq. (|127|l . is 



{nk')n 



Uk 

i n k} n() 



Pn ,{n k }n Q 
Pn 



(129) 



Therefore, the equation of motion for p no can be rewritten in the symmetrical and transparent form 



where 



jj:Pno = -n{K no (n + l)p no - K no -in Pn -i 
+H no n p no - H no+1 (n + l)p no+1 }, 



K n a = fa*' + 1 )( n fc')no- H no = ^' {{n k ' 



1). 



(130) 



(131) 



fc'>0 



fc'>0 



We can obtain the steady state distribution of the number of atoms condensed in the ground level of the trap from 
Eq. I|13(J|) • The mean value and the variance of the number of condensed atoms can then be determined. It is clear 
from Eq. (|130f) that there are two processes, cooling and heating. The cooling process is represented by the first two 
terms with the cooling coefficient K no while the heating by the third and fourth terms with the heating coefficient 
H no . The detailed balance condition yields the following expression for the number distribution of the condensed 
atoms 



Pn 



n 

Poll 



(132) 



where the partition function 



Zn — 



N N 

— =e n 



no— i— rio + 1 



Ki- 



(133) 



is determined from the normalization condition X)^ =oP™o = !• The functions Hi and Ki as given by Eq. (|131fl . 
involve, along with r\k> (Eq. H125(l ). the function (nk>)n (Eq. (|129Jl ). In the following sections, we shall derive closed 
form expressions for these quantities under various approximations. The master equation l|13(J|) for the distribution 
function for the condensed atoms is one of our main results. It yields explicit expressions for the statistics of the 
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FIG. 11: Physical interpretation of various coefficients in the master equations 



condensed atoms and the canonical partition function. Physical interpretation of various coefficients in the master 
equations is summarized in Fig. 1111 

Under the above assumption of a thermal equilibrium for non-condensed atoms, we have 



Mno = _L_lILo / ■ (134) 

In the next two sections we present different approximations that clarify the general result (|132ll . 

Short summary of this subsection is as follows. We introduce the probability of having n atoms in the ground level 
and n k atoms in the fc-th level P„ 0inil . »»*,...• We assume that the atoms in the excited levels with a given number of 
condensed atoms no arc in equilibrium state at the temperature T, then 

Pn , ni ,....n k .... = exp[-(3(E Q n + Ewi + ... + E k n k +)]. (135) 
An 



This equation yields 



P N ee P no=N . ni=0 ...., nk=0 .... = -^exp[-f3E N}. (136) 
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Assuming Eq = we obtain the following expression for the partition function 

Zn = (137) 
"n 

We assume that Bose gas can exchange heat (but not particles) with a harmonic-oscillator thermal reservoir. The 
reservoir has a dense and smooth spectrum. The average occupation number of the heat bath oscillator at a frequency 
u) q = cq is 

^ = — m 71 r\ T (138) 
ex-p(hu! q / kb-L ) — 1 

The master equation for the distribution function of the condensed bosons p na = p no .n takes the form 

Pno = -/] Kkq(n,k)n n {Vq + l )[( n Q + l)P«o ~ "-OPno-lkool- 
k.q 



- ^ K kq( n k + tynoVqlnoPna ~ (™0 + l)Pn + l]Hcat- (139) 

k,q 

The factors Kk q embody the spectral density of the bath and the coupling strength of the bath oscillators to the gas 
particles, and determine the rate of the condensate evolution since there is no direct interaction between the particles 
of an ideal Bose gas. Since Kk q — K ■ 5(Mlk — hcq) the sum J^k q reduces to 

-Pno = -y^(nk)n {Vk + l)[(»lQ + l)Pno ~ ^0Pn -l]Kool- 



'in Vk[n p na - (n + l)Pno+l]Heat- (140) 

k 

Particle number constraint comes in a simple way: X)fc( n fc)«o = N — fio- 



D. Low temperature approximation 

At low enough temperatures, the average occupations in the reservoir are small and rjk + 1 — 1 in Eq. 11311) . This 
suggests the simplest approximation for the cooling coefficient 

K no ~J2{ n k)n =N-n . (141) 

fe 

In addition, at very low temperatures the number of non-condensed atoms is also very small, we can therefore 
approximate (nk') no + 1 by 1 in Eq. (|131f) . Then the heating coefficient is a constant equal to the total average 
number of thermal excitations in the reservoir at all energies corresponding to the energy levels of the trap, 

H no ^H, W = £ffc=5>' h */ r -l)- 1 . (142) 

fc>0 fe>0 

Under these approximations, the condensate master equation l|13U|) simplifies considerably and contains only one 
non-trivial parameter Tl. We obtain 

^■Pn = -« {(A - n a )(n + l)p no - (A - n + l)n p no -i 

+U[n p no - (n + l)p„ 0+ i]}. (143) 

It may be noted that Eq. I|143|) has the same form as the equation (|1U7|) of motion for the photon distribution function 
in a laser operating not too far above threshold. The identification is complete if we define the gain, saturation, and 
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loss parameters in laser master equation by k(N + 1), k, and kTI, respectively. The mechanism for gain, saturation, 
and loss are however different in the present case. 

A laser phase transition analogy exists via the P-representation [Tfl Flj . The steady-state solution of the Fokker- 
Planck equation for laser near threshold is 0] 

P(a,a*) = lexp[(^)H 2 - ^|a| 4 ] (144) 

which clearly indicates a formal similarity between 

In P(a, a*) = - In Af + (1 - H/(N + l))n - (1/2(JV + l))n 2 (145) 

for the laser equation and the Ginzburg-Landau type free energy [T^, Fol Fl) 

G(no) = lnp no « const + a(T)no + b(T)n , (146) 

where \a\ 2 = n , a(T) = -(N - H)/N and b(T) = 1/(2N) for large N near T c . 

The resulting steady state distribution for the number of condensed atoms is given by 

I f^N—no 

Pn = -^-ITj \ii (!47) 

Z N (N - no)! 

where Zjv = 1/piv is the partition function. It follows from the normalization condition ^2 no p no = 1 that 

Z N = e n T(N + 1,H)/N\, (148) 

where T(a,x) — t a ~ 1 e~ t dt is an incomplete gamma-function. 

The distribution (|147|) can be presented as a probability distribution for the total number of non-condensed atoms, 
n = N — no , 

p ^ Pn -^yWTVh)^- (149) 

It looks somewhat like a Poisson distribution, however, due to the additional normalization factor, N\/T(N~\-1 1 Ti) ^ 1, 
and a finite number of admissible values of n = 0, 1, . . . , N, it is not Poissonian. The mean value and the variance 
can be calculated from the distribution l|147|) for an arbitrary finite number of atoms in the Bose gas, 

(n ) = N - H + H N+1 /Z N N\, (150) 



An% = <«o) - <«o) 2 = H (1 - ((n ) + 1)H N /Z N M) . (151) 

As we shall see from the extended treatment in the next section, the approximations (|141|> . Ijl42|l and, therefore, 
the results <|15U|) . i|151f) are clearly valid at low temperatures, i.e., in the weak trap limit, T <C E\, where E\ is an 
energy gap between the first excited and the ground levels of a single-particle spectrum in the trap. However, in the 
case of a harmonic trap the results (|150|) . (|151|) show qualitatively correct behavior for all temperatures, including 
T> £l and T-T c Hi. 

In particular, for a harmonic trap we have from Eq. I|142|l that the heating rate is 

" - 1» - £ exp^g + ' m + „) - 1] » ( 3 « 3 > - N (I) 3 ■ < 152 » 

Thus, in the low temperature region the master equation (|143fl for the condensate in the harmonic trap becomes 

(n + l) 2 ] p no + [(N + l)n - n 2 ]^-! 

-K + lK 0+ i]. (153) 



-Pn = - [(N + l)(n + l)- 
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E. Quasithermal approximation for non-condensate occupations 



At arbitrary temperatures, a very reasonable approximation for the average non-condensate occupation numbers 
in the cooling and heating coefficients in Eq. I|131|) is suggested by Eq. (|134fl in a quasithermal form, 



(n k ) no = Vk^2(n k ) no/ X] Vk> = ,^ /T Z^vu ' ( 154 ) 
fe>0 1 ' ' 1 



where e k = hv k , r] k is given by Eq. (|125fl and Ti by Eq. (|142|) . Equation 1|154|) satisfies the canonical-ensemble 
constraint, N = no + X)fc>o independently of the resulting distribution p no . This important property is based on 
the fact that a quasithermal distribution Ijl54(l provides the same relative average occupations in excited levels of the 
trap as in the thermal reservoir, Eq. (|125|l . 

To arrive at the quasithermal approximation in Eq. I|154fl one can g° along the following logic. In the low 
temperature limit we assumed rj k 1 and took 

^2(nk)n (Vk + 1) ~ ^2(n k )no =N-n . 

k k 

To go further, still in the low temperature limit, we can write 

exp(-(3E k ) 



(rik)n » (N - n ) 



Y,k>*M-PEk>) 



This is physically motivated since the thermal factor in [...] is the fraction of the excited atoms in the state k, and 
N ~ is the total number of excited atoms. Note that 



exp(/3£ fc )-l ^ rv ^ l + m' 

Since we are at low temperature we take exp(—(3E k ) ~ r\k and therefore 



/ \ ~ (at - \ Vk (N-no) 

(n k ) na ^(N-n )— = T — — — -, (155) 

J2 k Vk [exp(PE k ) - 1 ]H 



where 



k 

Now this ansatz is good for arbitrary temperatures. As a result, 

J2(n k ) no (Vk + l) » (N-n )(l + n), (156) 
k 

where 

Another line of thought is the following: 

rife 



(rik)no ~ (N - n )-. 



But by detailed balance in the steady state 

k(u + l)n k (fj k + 1) PS Kn (n k + l)fj k 

and if the ground level is macroscopically occupied then no ps no ± 1. Since even at T — T c one finds no ~ y/N, one 
"always" has fio 3> 1. Therefore, n k (fj k + 1) ~ (n k + l)fj k and, hence, n k pa fj k . As a result we again obtain Eq. (|155fl . 
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Calculation of the heating and cooling rates in this approximation is very simple. For example, for the heating rate 
we have 

J2^ n k + l )no(Vk) =£(%)+X>k)»o(%> ~n + V (N-n Q ). (157) 

In summary, the cooling and heating coefficients l|131(l in the quasithermal approximation of Eq. I|154|) are 

K no = (N-n ){l + ri), H no = H + (N - n )n. (158) 

Compared with the low temperature approximation 1141|) and i|142H . these coefficients acquire an additional contri- 
bution (N — no) 7 ? due to the cross-excitation parameter 

" = E^X^o = ^ E {e ej_ ir ( 159 ) 

u fe>0 k>0 v ; 

F. Solution of the condensate master equation 

Now, at arbitrary temperatures, the condensate master equation H13()(l contains two non-trivial parameters, Ji and 



1, 

" l>, '° = -k{(1 + r))[(N - n )(n + l) Pno -(N- n + l)n oPno _i 



dt 

+[H + (N — n Q )n]n oPno -[H + (N-n - l)r)]{n + l)p„ +i}- (160) 



It can be rewritten also in the equivalent form 
1 dp nQ 



= -[(N+ l)(n Q + 1) - (n + l) 2 }Pn + [(JV + l)n - n > no _i 
k at 



{T/T c f N[n oPno - (n + l>» 0+ i]. (161) 



The steady-state solution of Eq. I|160l) is given by 

\l + r>) Z N \ N-m / VI + 77/ ' { ' 



1 (N -n + H/ri-l)\ ( r\ \ N ~no _ l (N - n + % - V\ ( v ^-n 

Pn 



Z N {n/r]-iy.{N -n )l \1 + T)J Z N \ N-n /Vl + 77 

where the canonical partition function Z^ = 1/pat is 



N - n + H/v - 1\ / i] \ N ~ n ° 



N 

Ei n - n j vr 

n =0 V 



(163) 



It is worth noting that the explicit formula (|162fl satisfies exactly the general relation between the probability distri- 
bution of the number of atoms in the ground state, p no , and the canonical partition function |46j . Eq. I|79() . 

The master equation (jluOfl for p„ , and the analytic approximate expressions (|162f) and i|163|) for the condensate 
distribution function p no and the partition function Zn, respectively, are among the main results of the condensate 
master equation approach. As we shall see later, they provide a very accurate description of the Bose gas for a large 
range of parameters and for different trap potentials. Now we are able to present the key picture of the theory of 
BEC fluctuations, that is the probability distribution p no , Fig. 1121 Analogy with the evolution of the photon number 
distribution in a laser mode (from thermal to coherent, lasing) is obvious from a comparison of Fig. 1121 and Fig. 
EJ With an increase of the number of atoms in the trap, N, the picture of the ground-state occupation distribution 
remains qualitatively the same, just a relative width of all peaks becomes more narrow. 

The canonical partition function (|163|l allows us to calculate also the microcanonical partition function £l(E, N) by 
means of the inversion of the definition in Eq. (|74|) . Moreover, in principle, the knowledge of the canonical partition 
function allows us to calculate all thermodynamic and statistical equilibrium properties of the system in the standard 
way (see, e.g., and discussion in the Introduction). 

Previously, a closed-form expression for the canonical partition function (|74|l was known only for one-dimensional 
harmonic traps 68, 69] 



-m t ) = n x _^it - ( i64 ) 

k=l 



38 




FIG. 12: Probability distribution of the ground-state occupation, p„ , at the temperature T — 0.2T C in an isotropic harmonic 
trap with N = 200 atoms as calculated from the solution of the condensate master equation 11300 in the quasithermal 
approximation, Eq. (11621 . (solid line) and from the exact recursion relations in Eqs. 17911 and 1801 (dots). 



In the general case, there exists only the recursion relation H8U|) that is quite complicated, and difficult for analysis 

The distribution (|162[1 can also be presented as a probability distribution for the total number of non-condensed 
atoms, n — N — n , 

1 (n + H/r)-l\( V y 

Pn = PN-n = -S- ( T— • (165) 

Z N \ n J \ l + r)J 

The distribution (|165fl can be named as a finite negative binomial distribution, since it has the form of the well-known 
negative binomial distribution 

P„=(^ + ™~ 1 ^q n (l-q) M , n = 0,1, 2,..., oo, (166) 
that was so named due to a coincidence of the probabilities P n with the terms in the negative-power binomial formula 



1 ^ (n + M-l\ „ , 1cf7 , 

It has a similar semantic origin as the well-known binomial distribution, P n = ( )(1 — q) n q AI ~ n , which was named 

after a Newton's binomial formula [q+ (1 — q)] M = J2^=o ( M ) 0- ~ q) n q M ~ n - The finite negative binomial distribution 
(|165fl tends to the well-known distribution 1|166|) only in the limit N ^> (1 + ri)'H- 
The average number of atoms condensed in the ground state of the trap is 

N 

(no) = ^2 n p no . (168) 

"0=0 

It follows, on substituting for p no from Eq. (|162f> . that 

(n )=N-n+ Po r](N + n/r,) . (169) 

The central moments of the mth order, m > 1, of the number-of-condensed-atom and number-of- non-condensed- 
atom fluctuations are equal to each other for even orders and have opposite signs for odd orders, 

((no-ho) m ) = (-l) m ((n-nr). (170) 
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The squared variance can be represented as 

N 



An 2 , = (n 2 ) - (n) 2 = £ n(n - l)P n + (n) - (n) 2 (171) 



n=0 

and calculated analytically. We obtain 



An 2 



(1 + r])H - p (r)N + H) (N - H + 1 + 77) - pI^N + Hf . (172) 



where 



1 (N + H/r)-l)\ 



N 



Po = m{Hh-i)\ (173) 

is the probability that there are no atoms in the condensate. 

rr-^* N 

All higher central moments of the distribution Eq. i|162|) can be calculated analytically using (no s ) — n aP«o 

no=0 

and Eqs. (|162[1 . (|163fl . In particular, the third central moment is 

((n -(n » 3 ) = -(l + r ? )(l + 27 7 )W 

+ Po (H + t]N)[l + (H- N) 2 + 2(rf + N(l + rj)) + 3(v - H(l + 77))] 

+3pl(H + r]N) 2 (l +rj-H + N) + 2pl(H + f]N) 3 . (174) 

The first four central moments for the Bose gas in a harmonic trap with N = 200 atoms are presented in Fig. 1131 as 
the functions of temperature in different approximations. 

For the "condensed phase" in the thermodynamic limit, the probability po vanishes exponentially if the temperature 
is not very close to the critical temperature. In this case only the first term in Eq. I|172|) remains, resulting in 

An 2 = {l + V )H = Y,(( n k) 2 + (n k ))- (175) 

k>0 

This result was obtained earlier by standard statistical methods (see ^3| and references therein). 

It is easy to see that the result (|165ll reduces to the simple approximation l|149|) in the formal limit 77 — ^ 0, H./t] — > 
00, when 

W + H/v) ~* WJ ' (1?6) 

The limit applies to only very low temperatures, T e±. However, due to Eqs. 1142|) and i|159fl . the parameter Ti/rj 
tends to 1 as T — > 0, but never to infinity. Nevertheless, the results (|169fl and (117211 agree with the low temperature 
approximation results (|150fl and (|151fl for T -C £\. In this case the variance yAn^ is determined mainly by a square 
root of the mean value (n) which is correctly approximated by Eq. (|150|l as (n) = N — (n ) » H. 



G. Results for BEC statistics in different traps 



As we have seen, the condensate fluctuations are governed mainly by two parameters, the number of thermal 
excitations TL and the cross-excitation parameter rj. They are determined by a single-particle energy spectrum of the 
trap. We explicitly present them below for arbitrary power-law trap. We discuss mainly the three-dimensional case. 
A generalization to other dimensions is straightforward and is given in the end of this subsection. First, we discuss 
briefly the case of the ideal Bose gas in a harmonic trap. It is the simplest case since the quadratic energy spectrum 
implies an absence of the infrared singularity in the variance of the BEC fluctuations. However, because of the same 
reason it is not robust relative to an introduction of a realistic weak interaction in the Bose gas as is discussed in 
Section V. 
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FIG. 13: The first four central moments for the ideal Bose gas in an isotropic harmonic trap with N = 200 atoms as calculated 
via the solution of the condensate master equation (solid lines - quasithermal approximation, Eq. (11621 : dashed lines - low 
temperature approximation, Eq. an d via the exact recursion relations in Eqs. I|79|l and ilSIH (dots). 



Harmonic trap 



The potential in the harmonic trap has, in general, an asymmetrical profile in space, V ext (x, y, z) = ^(x 2 u> 2 
y 2 ijJy + z 2 w 2 ), with eigenfrequencies {uj x ,ujy,uj z } = u>, uj x > aj y > to x > 0. Here m is the mass of the atom. The 
single-particle energy spectrum of the trap, 

£k = hhuj = h(k x uj x + kyUiy + k z u> z ), (177) 

can be enumerated by three non-negative integers {k x , k y , k z } = k, k XjV ^ z > 0. We have 

H = 12 jaJ/T _ 1 ■ ^ H = 12 ( e Kk»/T _ D2 ■ ( 178 ) 

k>0 k>0 ^ ' 

The energy gap between the ground state and the first excited state in the trap is equal to e\ = frus x . 

If the sums can be replaced by the integrals (continuum approximation), i.e., if ftw x <C T, the parameters Ti. and 
r/Ti. are equal to 

T 3 / T \ 3 

K= T3 C(3 = (=-) N, (179) 



rfli 



rp3 



(C(2)-C(3)) = © 



H 3 UJ x LOyUJ z 



C(3) 



C(3) 



(180) 
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where a standard critical temperature is introduced as 

T ^ h (^Wr) 1/3 -' C0) = 1.202..., C (2) = £. (181) 

Therefore, the cross-excitation parameter rj is a constant, independent of the temperature and the number of atoms, 
j] = [C(2) - C(3)]/C(3) ~ 0.37. The ratio H/r] = iV(T/T c ) 3 [C(3)/(C(2) - C(3))] goes to infinity in the thermodynamic 
limit proportionally to the number of atoms TV. 

In the opposite case of very low temperatures, T <C huj x , we have 

H«exp(-^)+exp(-^)+exp(-^), (182) 

27^ 2^ 2/iw 2 
7 ? Wwcxp( — ) + exp( ^)+exp( — ) (183) 

with an exponentially good accuracy. Now the cross-excitation parameter 77 depends exponentially on the temperature 
and, instead of the number 0.37, is exponentially small. The ratio 

H = [exp(-^) + exp(-^) + exp(-^)] 2 ^ i 
V cxp(-2^) + exp(-^) + cxp(-^) ~ 

becomes approximately a constant. The particular case of an isotropic harmonic trap is described by the same 
equations if we substitute u> x = uj v — u> z = to. 

2. Arbitrary power-law trap 

We now consider the general case of a d-dimensional trap with an arbitrary power-law single-particle energy spectrum 

d 

e k = tiJ2^J, k {/.-,: j = l,2,...,d}, (185) 
j=i 

where kj > is a non-negative integer and a > is an index of the energy spectrum. We assume < w\ < L02 < 
• ■ • < ^di so that the energy gap between the ground state and the first excited state in the trap is E\ = hw\. We then 
have 

U = Y.^F— V ^ = E (e e k /T_ 1)2 - (186) 

k>0 k>0 v ' 

In the case E\ <§; T, the sum can be replaced by the integral only for the parameter TC (Eq. (|186fl ) if d > a, 

H = A^)T d /° = (y)^^ d>a > (187) 

where the critical temperature is 



N 



*/* A n- n + 1)] 



d 1. , \ l / a 



.^C(d/«r). 

The second parameter can be calculated by means of this continuum approximation only if < a < d/2, 



V H = AT d l a (C(^ - 1) - C(|)) = (§) C( " } , < a < d/2. (189) 



If (j > d/2, it has a formal infrared divergence and should be calculated via a discrete sum, 

T \ 2 ^r2ald a <y-d 



r,H=(±) N 2 °/ d ^ z-n, a>d/2, (190) 

\tJ [r(i + i)] 2CT [ C (|)] 2CT/<i 
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where 



a a ,d - 



k>0 



k 



The traps with the dimension lower than the critical value, d < a, can be analyzed on the basis of Eqs. (|186fl as well. 
We omit this analysis here since there is no phase transition in this case. 

The cross-excitation parameter r\ has different dependence on the number of atoms for high, d > 2a, or low, d < 2a, 
dimensions, 

C(d _ i)- ((d) 

71 C (i) ' d >^>^ (191) 



" ' N^— ,S — „ , d<2a. (192) 



[r(i + i)] ZCT [C(|)]- 



Therefore, the traps with small index of the energy spectrum, < a < d/2, are similar to the harmonic trap. The 
traps with larger index of the energy spectrum, a > d/2, are similar to the box with "homogeneous" Bose gas. For the 
latter traps, the cross-excitation parameter n goes to infinity in the thermodynamic limit, proportionally to N 2a / d ^ 1 . 
The ratio TL/rj goes to infinity in the thermodynamic limit only for < a < d. In the opposite case, a > d, it goes to 
zero. We obtain 

- = (t) N 7Td i\ 77*-V d>2a>0, (193) 



- = (-) N ^-°M r - + 1 c - 

It is remarkable that BEC occurs only for those spatial dimensions, d > a, for which 7i/r) — * oo at iV — * oo. (We 
do not consider here the case of the critical dimension d = a, e.g., one-dimensional harmonic trap, where a quasi- 
condensation occurs at a temperature T c ~ huiN/ log AT.) For spatial dimensions lower than the critical value, d < a, 
BEC does not occur (see, e.g., ). Interestingly, even for the latter case there still exists a well-defined single peak 
in the probability distribution p no at low enough temperatures. With the help of the explicit formulas in Section III 
we can describe this effect as well. 

In the opposite case of very low temperatures, T <C £i, the parameters 

d d 2hu, 

H«^e^, ?77i«^e - ^", rj - (195) 

i=i j=i 

are exponentially small. The ratio ^ ~ Sj=i e _ ^ riWj_£1 ^ T ~ d becomes a constant. 

Formulas p85(l - (|195|l for the arbitrary power-law trap contain all particular formulas for the three-dimensional 
harmonic trap (d = 3, a = 1) and the box, i.e., the "homogeneous gas" with d = 3 and a = 2, as the particular 
cases. 

In Fig. EH numerical comparison of the results obtained from the exact recursion relations in Eqs. H79 M 80 |I an< l 
from our approximate explicit formulas from Section IV in the particular case of the ideal Bose gas in the three- 
dimensional isotropic harmonic trap for various temperatures is demonstrated. The results indicate an excellent 
agreement between the exact results and the results based on quasithermal approximation, including the mean value 
(no)) the squared variance Anp as well as the third and fourth central moments. The low temperature approximation, 
Eq. (|147fl . is good only at low temperatures. That is expected since it neglects by the cross-excitation parameter rj. 



4(7 1 a 



-1 
a.d 



d < 2a. 



(194) 



H. Condensate statistics in the thermodynamic limit 

The thermodynamic, or bulk limit implies an infinitely large number of atoms, — > oo, in an infinitely large 
trap under the condition of a fixed critical temperature, i.e., Nu) x oj y uj z = const in the harmonic trap, L 3 N = const 
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in the box, and N°"IIj_ 1 u)j = const in an arbitrary d-dimensional power-law trap with an energy spectrum index 
a. Then, BEC takes place at the critical temperature T c (for d > a) as a phase transition, and for some lower 
temperatures the factor po is negligible. As a result, we have the following mean value and the variance for the 
number of condensed atoms 

( no )=N-H = N- y £ l ^- l , (196) 

fe>0 



Ang = (1 + rj)H ^ £ + ^ , e J_ 1)2 , (197) 

fe>0 fc>0 v ; 

which a gree with the results obtained for the ideal Bose gas for different traps in the canonical ensemble by other 
authors jl3t [Tj, |46|, [73|, |7J, |75|, IZa |Z3 ■ I* 1 particular, we find the following scaling of the fluctuations of the number 
of condensed atoms: 

/ \ d/a 
A^r) N , d>2o">0x 

Ang~Cx V c/ 2 , ei «T<r c , (198) 

V (f) N^/ d , d < 2a / 



Ang 



exp{- 



n 



1 1 / d 



a/d 



1) 



TN a / d 



}, T« £l , 



(199) 



where C is a constant. From Eq. (|198|) . we see that in the high dimensional traps, d > 2a, e.g., in the three-dimensional 
harmonic trap, fluctuations display the proper thermodynamic behavior, Atiq cx N. However, fluctuations become 
anomalously large 0, EH 0> El] > An§ cx N 2cr / d » N, in the low dimensional traps, a < d < 2a. In the quantum 
regime, when the temperature is less than the energy gap between the ground and the first excited level in the trap, 
it follows from Eq. I|199(l that condensate fluctuations become exponentially small. For all temperatures, when BEC 
exists (d > a) , the root-mean-square fluctuations normalized to the mean number of condensed atoms vanishes in the 
thermodynamic limit: yf Ang/ (no) — > as N — > oo. 

Another remarkable property of the distribution function obtained in Section IV is that it yields the proper mean 
value and variance of the number of atoms in the ground level of the trap even for temperatures higher than the 
critical temperature. In particular, it can be shown that its asymptotic for high temperatures, T S> T c , yields the 
standard thermodynamic relation Ano ~ (hq) known from the analysis of the grand canonical ensemble [13j- This 
nice fact indicates that the present master equation approach to the statistics of the cool Bose gas is valuable in the 
study of mesoscopic effects as well, both at T < T c and T > T c . Note that, in contrast, the validity of the Maxwell's 
demon ensemble approach j^j to the statistics of BEC remains restricted to temperatures well below the onset of 
BEC, T <T C . 



I. Mesoscopic and dynamical effects in BEC 

In recent experiments on BEC in ultracold gases pil l25l l26l l27t I2H , |29| , the number of condensed atoms in the trap 
is finite, i.e., mesoscopic rather than macroscopic, N ~ 10 3 — 10 6 . Therefore, it is interesting to analyze mesoscopic 
effects associated with the BEC statistics. 

The mean number of atoms in the ground state of the trap with a finite number of atoms is always finite, even 
at high temperatures. However, it becomes macroscopically large only at temperatures lower than some critical 
temperature, T c , that can be defined via the standard relation 

J2vk(T c ) = H(T c ) = N. (200) 

This equation has an elementary physical meaning, namely it determines the temperature at which the total average 
number of thermal excitations at all energy levels of the trap becomes equal to the total number of atoms in the 
trap. The results (|162fl . I|lt)9l) . I|172|) shown in Fig. 1131 explicitly describe a smooth transition from a mesoscopic 
regime (finite number of atoms in the trap, N < 00) to the thermodynamic limit (N — 00) when the threshold of 
the BEC becomes very sharp so that we have a phase transition to the Bose-Einstein condensed state at the critical 



44 



temperature given by Eq. H200|) . This can be viewed as a specific demonstration of the commonl y ac cepted resolution 
to the Uhlenbeck dilemma in his famous criticism of Einstein's pioneering papers on BEC [a. l9l floL fl^ | . 

Although for systems containing a finite number of atoms there is no sharp critical point, as is obvious from Figs. 03 
1121 and El it is useful to define a critical characteristic value of a temperature in such a case as well. It should 
coincide with the standard definition 12001) in the thermodynamic limit. Different definitions for T c were proposed 
and discussed in 0, 0, H3, El l82t l83l Is3 . Isfij . We follow a hint from laser physics. There we know that 
fluctuations dominate near threshold. However, we define a threshold inversion as that for which gain (in photon 
number for the lasing mode) equals loss. Let us use a similar dynamical approach for BEC on the basis of the master 
equation, see also |86j . 

We note that, for a laser operating near the threshold where B/A -C 1, the equation 1118fl of motion for the 
probability p n of having n photons in the cavity implies the following rate of the change for the average photon 
number: 

^(n) = (A-C)(n)-B((n + l) 2 )+A. (201) 

Here A, B, and C are the linear gain, nonlinear saturation, and linear loss coefficients, respectively. On neglecting the 
spontaneous emission term A and noting that the saturation term B((n + l) 2 ) is small compared to (A — C)(n) near 
threshold, we define the threshold (critical) inversion to occur when the linear gain rate equals the linear loss rate, 
i.e., A = C. 

Similar to laser physics, the condensate master equation l|130|) implies a coupled hierarchy of moment equations 
which are useful in the analysis of time evolution. In the quasithermal approximation i|160|) , we find 

d ( n o) _ f l \sn j_ m \mf/*,i\ _l _ _ /^+2\ 



dt 



K E0 {(1+ ^ ((n ° ) + (n ° +1)) - (R 



/ \"0 



+(-iy-\n + vn)H +x ) - (-i) l - i v {4 +2 )}. (202) 

Similar moment equations in the low-temperature approximation (|143fl follow from Eq. I|202|) with r\ = Q, 



d(n l ) J^fl 



dt 

i=0 



«E (|W<»o> + W ~ 1)K +1 > ~ K +2 > + (-l)'-*?^ 1 )}. (203) 

i— n \ / 



The dynamical equation for the first moment, as follows from Eq. I|202|) . has the following form 

d(n ) 



dt 



= «{(1 + V )N + (N - 1 - n - H)(n ) - (n 2 )}. (204) 



Near the critical temperature, T « T c , the mean number of the condensed atoms is small, (no) <C N, and it is 
reasonable to neglect the second moment (n 2 ,) compared to N(n ) and the spontaneous cooling (spontaneous emission 
in lasers) term k(1 + 77)^ compared to kJV(tiq). In this way, neglecting fluctuations, we arrive at a simple equation 
for the competition between cooling and heating processes, 

-M««(iV_W-q)(no). (205) 

In analogy with the laser threshold we can define the critical temperature, T = T c , as a point where cooling equals 
heating, i.e., d(no)/dt = 0. This definition of the critical temperature 

H(T C ) + V {T C ) = N, (206) 

is valid even for mesoscopic systems and states that at T = T c the rate of the removal of atoms from the ground 
state equals to the rate of the addition, in the approximation neglecting fluctuations. In the thermodynamic limit 
it corresponds to the standard definition, Eqs. (|181(l and l|188[) . For a mesoscopic system, e.g., of N = 10 3 atoms 
in a trap, the critical temperature as given by Eq. (|2U6|I is about few per cent shifted from the thermodynamic- 
limit value given by Eqs. 118111 and i|188l) . Other definitions for T c also describe the effect of an effective-T c shift 
[TTL l66l Ff^. l&oL l8ll 18^. l83l l84l l85| . which is clearly seen in Fig. 1141 and agree qualitatively with our definition. 

Note that precisely the same definition of the critical temperature follows from a statistical mechanics point of view, 
which in some sense is alternative to the dynamical one. We may define the critical temperature as the temperature 
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at which the mean number of condensed atoms in the steady- state solution to the master equation vanishes when 
neglecting fluctuations and spontaneous cooling. We make the replacement (rig) ~ ( n o) 2 m Eq. I|2U4|) and obtain the 
steady-state solution to this nonlinear equation, (no) = N — W — 77. Now we see that (no) vanishes at the same critical 
temperature (|2U6[) . Finally, we remind again that a precise definition of the critical temperature is not so important 
and meaningful for the mesoscopic systems as it is for the macroscopic systems in the thermodynamic limit since for 
the mesoscopic systems, of course, there is not any sharp phase transition and an onset of BEC is dispersed over a 
whole finite range of temperatures around whatever T c , as is clearly seen in Figs. |3 El and 1 141 

V. QUASIPARTICLE APPROACH AND MAXWELL'S DEMON ENSEMBLE 

In order to understand relations between various approximate schemes, we formulate a systematic analysis of the 
equilibrium canonical-ensemble fluctuations of the Bosc-Einstein condensate based on the particle number conserv- 
ing operator formalism of Girardeau and Arnowitt |8?|. and the concept of the canonical-ensemble quasiparticles 
|2Ct l2l| . The Girardeau- Arnowitt operators can be interpreted as the creation and annihilation operators of the 
canonical-ensemble quasiparticles which are essentially different from the standard quasiparticles in the grand canon- 
ical ensemble. This is so because these operators create and annihilate particles in the properly reduced many-body 
Fock subspace. In this way, we satisfy the TV-particle constraint of the canonical-ensemble problem in Eq. I|72[) from 
the very beginning. Furthermore, we do this while taking into account all possible correlations in the TV-boson system 
in addition to what one has in the grand canonical ensemble. These canonical-ensemble quasiparticles fluctuate inde- 
pendently in the ideal Bose gas and form dressed canonical-ensemble quasiparticles in the dilute weakly interacting 
Bose gas due to Bogoliubov coupling (see Section VI below). 

Such an analysis was elaborated in [23, and resulted in the explicit expressions for the characteristic function 
and all cumulants of the ground-state occupation statistics both for the dilute weakly interacting and ideal Bose gases. 
We present it here, including the analytical formulas for the moments of the ground-state occupation fluctuations 
in the ideal Bose gases in an arbitrary power-law trap, and, in particular, in a box ("homogeneous gas") and in an 
arbitrary harmonic trap. In Section VI we extend this analysis to the interacting Bose gas. In particular, we calculate 
the effect of Bogoliubov coupling between quasiparticles on suppression of the ground-state occupation fluctuations 
at moderate temperatures and their enhancement at very low temperatures and clarify a crossover between ideal-gas 
and weakly-interacting-gas statistics which is governed by a pair-correlation, squeezing mechanism. The important 
conclusion is that in most cases the ground-state occupation fluctuations are anomalously large and are not Gaussian 
even in the thermodynamic limit. 

Previous studies focused mainly on the mean value, no, and squared variance, ((no — no) 2 ), of the number of 
condensed atoms [8^. Higher statistical moments are more difficult to calculate, and it was often assumed that 
the condensate fluctuations have vanishing higher cumulants (semi- invariants) . That is, it was assumed that the 
condensate fluctuations are essentially Gaussian with all central moments determined by the mean value and the 
variance. We here show that this is not true even in the thermodynamic limit. In, particular, we prove that in the 
general case the third and higher cumulants normalized by the corresponding power of the variance do not vanish 
even in the thermodynamic limit. 

The results of the canonical-ensemble quasiparticle approach are valid for temperatures a little lower than a critical 
temperature, namely, when the probability of having zero atoms in the ground state of the trap is negligibly small and 
the higher order effects of the interaction between quasiparticles are not important. We outline also the Maxwell's 
demon ensemble approximation introduced and studied for the ideal Bose gas in 0> 13 03 EH and show that 
it can be justified on the basis of the method of the canonical-ensemble quasiparticles, and for the case of the ideal 
Bose gas gives the same results. 

This Section is organized as follows: We start with the reduction of the Hilbert space and the introduction of the 
canonical-ensemble quasiparticles appropriate to the canonical-ensemble problem in subsection A. Then, in subsection 
B, we analytically calculate the characteristic function and all cumulants of the ground-state occupation distribution 
for the ideal Bose gas in a trap with an arbitrary single-particle energy spectrum. We also discuss the Maxwell's 
demon ensemble approach and compare it with the canonical ensemble quasiparticle approach. In subsection C we 
apply these results to the case of an arbitrary d-dimcnsional power-law trap which includes a three-dimensional box 
with periodic boundary conditions ( "homogeneous gas" ) and a three-dimensional asymmetric harmonic trap as the 
particular cases. 
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A. Canonical-ensemble quasiparticles in the reduced Hilbert space 

In principle, to study the condensate fluctuations, we have to fix only the external macroscopical and global, 
topological parameters of the system, like the number of particles, temperature, superfluid flow pattern ("domain" or 
vortex structure), boundary conditions, etc. We then proceed to find the condensate density matrix via a solution of 
the von Neumann equation with general initial conditions admitting all possible quantum states of the condensate. 
In particular, this is a natural way to approach the linewidth problem for the atom laser [89l l90| . Obviously, this is 
a complicated problem, especially for the interacting finite-temperature Bose gas; because of the need for an efficient 
technique to account for the additional correlations introduced by the constraint in such realistic ensembles. The 
latter is the origin of the difficulties in the theory of the canonical or microcanonical ensembles (see discussion in the 
Section III. A). According to 32], a calculation of equilibrium statistical properties using the grand canonical ensemble 
and a perturbation series will be impossible since the series will have zero radius of convergence. 

One way out of this problem is to develop a technique which would allow us to make calculations in the constrained 
many-body Hilbert space, e.g., on the basis of the master equation approach as discussed in Section IV. Another 
possibility is to solve for the constraint from the very beginning by a proper reduction of the many-body Hilbert 
space so that we can work with the new, already unconstrained quasiparticles. This approach is demonstrated in the 
present Section. Working in the canonical ensemble, we solve for the fluctuations of the number of atoms in the ground 
state in the ideal Bose gas in a trap (and similarly in the weakly interacting Bose gas with the Bogoliubov coupling 
between excited atoms, see Section VI). More difficult problems involving phase fluctuations of the condensate with 
an accurate account of the quasiparticle renormalization due to interaction at finite temperatures and the dynamics 
of BEC will be discussed elsewhere. 

We begin by defining an occupation number operator in the many-body Fock space as usual, 

n k = a+a k , u k |^ k n) )=n|V k n) ). a+\^ } ) =Vn~+T\^ +1) ). (207) 

The particle number constraint l|72|) determines a canonical-ensemble (CE) subspace of the Fock space. Again we 
would like to work with the particle-number conserving creation and annihilation operators. The latter are given in 
the Girardeau and Arnowitt paper |87j . 

/3k = /3 + a k , $o = (1 + n r 1/2 a . (208) 

These operators for k ^ can be interpreted as describing new canonical-ensemble quasiparticles which obey the 
Bose canonical commutation relations on the subspace no ^ 0, 

(209) 

We are interested in the properties of the fraction of atoms condensed in the ground level of the trap, k = 0. We focus 
on the important situation when the ground-state occupation distribution is relatively well peaked, i.e., its variance 
is much less than the mean occupation of the ground level of a trap, 

((n - n ) 2 ) 1/2 < n . (210) 

In such a case, the relative role of the states with zero ground-state occupation, hq = 0, is insignificant, so that we can 
approximate the canonical-ensemble subspace Tt by the subspace 7i^^L. Q . Obviously, this approximation is valid 
only for temperatures T < T c . 

The physical meaning of the canonical-ensemble quasiparticles, /3 k = /3o~a k , is that they describe transitions between 
ground (k = 0) and excited (k ^ 0) states. All quantum properties of the condensed atoms have to be expressed via 
the canonical-ensemble quasiparticle operators in Eq. I|208f) . In particular, we have the identity 

n =JV-^n k) (211) 

k^O 

where the occupation operators of the excited states are 

n k = a+a k = /3+/3 k . (212) 

Note that in Refs. quasiparticle operators similar in spirit to those of Ref. [87j were introduced which, unlike 

/3 k , did not obey the Bose commutation relations (|209|l exactly, if non-commutation of the ground-state occupation 
operators fig and is important. As was shown by Girardeau |93| . this is important because the commutation 
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corrections can accumulate in a perturbation series for quantities like an S'-matrix. Warning concerning a similar 
subtlety was stressed some time ago [94|. 

We are interested in fluctuations in the number of atoms condensed in the ground state of a trap, no- This is 
equal to the difference between the total number of atoms in a trap and the number of excited atoms, no = N — n. 
In principle, a condensed state can be defined via the bare trap states as their many-body mixture fixed by the 
interaction and external conditions. Hence, occupation statistics of the ground as well as excited states of a trap 
is a very informative feature of the BEC fluctuations. Of course, there are other quantities that characterize BEC 
fluctuations, e.g., occupations of collective, dressed or coherent excitations and different phases. 

B. Cumulants of BEC fluctuations in an ideal Bose gas 

Now we can use the reduced Hilbert space and the equilibrium canonical-ensemble density matrix p to conclude that 
the occupation numbers of the canonical-ensemble quasiparticles, n k , k ^ 0, are independent stochastic variables 
with the equilibrium distribution 

PkK) = exp(-n k e k /T)(l - exp(-e k /T)). (213) 

The statistical distribution of the number of excited atoms, n — X^k^o 7 ^' which is equal, according to Eq. (|211|l . 
to the number of non-condensed atoms, is a simple "mirror" image of the distribution of the number of condensed 
atoms, 

p(n) = p (n =N-n). (214) 
A useful way to find and to describe it is via the characteristic function 

Q n (u) =Tr{e iuA p}. (215) 
Thus upon taking the Fourier transform of 0„(it) we obtain the probability distribution 

P(n) = 7^ / e~ mn Q n (u)du. (216) 

Taylor expansions of & n (u) and logO„(u) give explicitly initial (non-central) moments and cumulants, or semi- 
invariants |72l l95| : 

e n (u) = a ™^> a ™ = < n ™> = 77^r e »( u )l«=0' ( 217 ) 

TTX . (XU 
log0„(u) = V K m ^-, K m = ——\oge n (u)\ u=0 , &n(u = 0) = 1. (218) 

— ' ml d{iu) m 

m— 1 v ' 

The cumulants K r , initial moments a m , and central moments p m = ((n — n) m ) are related to each other by the simple 
binomial formulas [7^. l95| via the mean number of the non-condensed atoms n — N — uq, 

Pr = ^(-lf JdrV, a r = Y^ (IW-fe^, 

fc=0 ^ ' k=0 ^ ' 

fi = Ki, ((71 - fif) = p 2 — K2, ((n ~ n) 3 ) = /x 3 = k 3 , ((n - n) 4 ) = pi — k 4 + 3^2, 



((n — n) 5 ) = p 5 = k 5 + 10k 2 k 3 , ((n — n) 6 ) = p 6 = k 6 + 15k 2 (k 4 + K 2 ) + 10k 3 , . . . (219) 

Instead of calculation of the central moments, p m = ((n — n) m ), it is more convenient, in particular so for the 
analysis of the non-Gaussian properties, to solve for the cumulants K m , which are related to the moments by simple 
binomial expressions. The first six are 

K\ =n, n 2 = p 2 , k 3 = p 3 , K 4 = p 4 - 3p%, 
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«5 = A*5 - 10M2M3, «6 = - 15/^2 (A*4 - 2jU§). (220) 

As discussed in detail below, the essence of the BEC fluctuations and the most simple formulas are given in terms of 
the "generating cumulants" k m which are related to the cumulants n m by the combinatorial formulas in Eq. 



Ki = ki, k 2 — k 2 + hi, k 3 = k 3 + 3/*2 + K 4 — ^4 + 6K3 + 7k,2 + ki, ■ ■ ■ (221) 

The main advantage of the cumulant analysis of the probability distribution p(n) is the simple fact that the 
cumulant of a sum of independent stochastic variables is equal to a sum of the partial cumulants, n r = X^k^o K r- 
This is a consequence of the equalities log0„(w) = logilk^o@n k (ii) = ^2^ ^og& nk (u) . For each canonical-ensemble 
quasiparticle, the characteristic function can be easily calculated from the equilibrium density matrix as follows 

Q nk (u) = Tr{e iuftk p k } = Tr{e ltlflk e~ £k " k/T }fl - e~ £k/T \ = fii-li (222) 

V / Zk — z 

Here we introduced the exponential function of the single-particle energy spectrum £k, namely Zk = ex p(£k/? 1 ), and 
a variable z = exp(iu) which has the character of a "fugacity" . As a result, we obtain an explicit formula for the 
characteristic function and all cumulants of the number of excited (and, according to the equation no — N — n, 
condensed) atoms in the ideal Bose gas in an arbitrary trap as follows: 

2k -1\ {e m -l) m ^ (iu) r 



log6 n (u) = ^log(^— -) = * 



.z\r — zJ ' — ' ml * — ' r! 

k^o m=i r=i 



(m - 1)! E( e£k/T - *r = E 4 m) «m- (223) 

k^o m=l 



Here we use the Stirling numbers of the 2nd kind |7^. 

- m / \ oo 

= h E(-D m - fe (T) - 1)" = « E 

t — n \ / „ i-. 



(224) 



that yield a simple expression for the cumulants n r via the generating cumulants k m . In particular, the first four 
cumulants are given in Eq. (|221|) . 

Thus, due to the standard relations H219|l . the result (|223Jl yields all moments of the condensate fluctuations. Except 
for the average value, all cumulants are independent of the total number of atoms in the trap; they depend only on the 
temperature and energy spectrum of the trap. This universal temperature dependence of the condensate fluctuations 
was observed and used in |42l [H . I4H IH4I l57j to study the condensate fluctuations in the ideal Bose gas on the basis 
of the so-called Maxwell's demon ensemble approximation. The method of the canonical-ensemble quasiparticles also 
agrees with and provides further justification for the "demon" approximation. The main point is that the statistics 
is determined by numbers and fluctuations of the excited, non-condensed atoms which behave independently of the 
total atom number iV for temperatures well below the critical temperature since all "excess" atoms stay in the ground 
state of the trap. Therefore, one can calculate statistics in a formal limit as if we have an infinite number of atoms in 
the condensate. That is we can say that the condensate plays the part of an infinite reservoir for the excited atoms, 
in agreement with previous works [Tj, H3, 0, Hy, |5J, |57j ■ 

Obviously, our approximation in Eq. I|21U|) as well as the Maxwell's demon ensemble approximation does not 
describe all mesoscopic effects that can be important very close to the critical temperature or for a very small number 
of atoms in the trap. However, it takes into account the effect of a finite size of a trap via the discreteness of the energy 
levels £k, i.e., in this sense the approximation (|210|) describes not only the thermodynamic limit but also systems with 
a finite number of atoms N. In addition, the mesoscopic effects can be partially taken into account by a "grand" 
canonical approximation for the probability distribution of the canonical-ensemble quasiparticle occupation numbers 

Pk(rik) = exp(-n k £k/7 1 )(l - exp(-e k /7 1 )), e k = £k - Mi (225) 

where the chemical potential is related to the mean number of the condensed atoms no = 1/(1 — exp(— (3p)) and 
should be found self-consistently from the grand-canonical equation 



k^O 



(226) 
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The canonical-ensemble quasiparticle result for all cumulants remains the same as is given by Eq. (|223|l above, with 
the only difference that now all quasiparticle energies are shifted by a negative chemical potential (e^ = — /i), 

r 

~ Km = (m-l)\^(e Sk / T -l)~ m , m = l,2,...; K r = 4^h m . (227) 

k^O m=l 

The first, m = 1, equation in Eq. I|227|) is a self-consistency equation l|22t)|) . The way it takes into account the 
mesoscopic effects (within this mean-number "grand" canonical approximation) is similar to the way in which the 
self-consistency equation (I264f) of the mean-field Popov approximation takes into account the effects of weak atomic 
interaction. The results of this canonical-ensemble quasiparticle approach within the "grand" canonical approximation 
for the quasiparticle occupations l|225l) were discussed in Section III for the case of the isotropic harmonic trap. 
Basically, the "grand" canonical approximation improves only the result for the mean number of condensed atoms 
no(T), but not for the fluctuations. 



C. Ideal gas BEC statistics in arbitrary power-law traps 

The explicit formulas for the cumulants demonstrate that the BEC fluctuations depend universally and only on the 
single-particle energy spectrum of the trap, £k- There are three main parameters that enter this dependence, namely, 

(a) the ratio of the energy gap between the ground level and the first excited level in the trap to the temperature, 

ei/T, 

(b) the exponent of the energy spectrum in the infrared limit, oc k a at k — > 0, and 

(c) the dimension of the trap, d. 

The result (|223|l allows to easily analyze the condensate fluctuations in a general case of a t rap with an arbitrary 
dimension d > 1 of the space and with an arbitrary power-law single-particle energy spectrum |4fii l55l l73| 

ei = hY / u j lj, l = {lf, j = l,2,...,d}, (228) 
i=i 

where lj > is a non-negative integer and a > is an index of the energy spectrum. The results for the particular 
traps with a trapping potential in the form of a box or harmonic potential well can be immediately deduced from the 
general case by setting the energy spectrum exponent to be equal to a — 2 for a box and a = 1 for a harmonic trap. 
We assume that the eigenfrequencies of the trap are ordered, < lu\ < loi < ■ ■ ■ < tod, so that the energy gap between 
the ground state and the first excited state in the trap is E\ — hw\. All cumulants i|223|) of the condensate occupation 
fluctuations are given by the following formula: 

f- d r 

K m = (m - 1)! Y, NtyE^ 1 ") - 1 ]~ m ' Kr = E fJ r m) ^' ( 229 ) 

\=(h,...J d )>0 3=1 m=l 

Let us consider first again the case of moderate temperatures larger than the energy gap, E\ ^ T < T c . The first 
cumulant, i.e., the mean number of non-condensed atoms, can be calculated by means of a continuum approximation 
of the discrete sum by an integral if the space dimension of a trap is higher than a critical value, d > a. Namely, one 
has an usual BEC phase transition with the mean value 

ki =n = N-n = A(^T d / a = (yY^ N ' d > ^ e x <^T<T c , (230) 
where the standard critical temperature is 



N 



A((d/a) 



a/d T(I + \)\ d 

A=^± 231 



The second-order generating cumulant can be calculated by means of this continuum approximation only if d > 2a, 
„ = AT'" U ( i - 1) - C ( <j) = (IY"N «1 - ? - «*> , (232) 



a a 



Tc> C(f) 
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In the opposite case it has to be calculated via a discrete sum because of a formal infrared divergence of the integral. 
Keeping only the main term in the expansion of the exponent in Eq. I|229|l . exp(y ^jlj) — 1 ~ T S-*=i we 

find 

K2 = {^' d /[T^ + l)r[a-)r' d ?a^ a<d<2a, (233) 
where — J2\>o(^ < j=i^ UJ j) 2 ^ d l £ i - The ratio of the variance to the mean number of non-condensed atoms is equal 



to v / K^/Ki = WK 1 + fta/ftf, i.e., 



^^^-y^'^K-m^, (234) 



a <d<2a. (235) 



We see that the traps with a relatively high dimension of the space, d > 2a, produce normal thermodynamic 
fluctuations (|234|l cx N^ 1 ' 2 and behave similar to the harmonic trap. However, the traps with a relatively low 
dimension of the space, a < d < 2a, produce anomalously large fluctuations (|235() in the thermodynamic limit, 
oc N /^ 1 3> N~ x / 2 and behave similar to the box with a homogeneous Bose gas, where there is a formal infrared 
divergence in the continuum-approximation integral for the variance. 

The third and higher-order central moments ((n — n) m ), or the third and higher-order cumulants n m , provide 
further parameters to distinguish different traps with respect to their fluctuation behavior. The m-th order generating 
cumulant can be calculated by means of the continuous approximation only if d > ma, 

— — dt = — — -, -r- / -7-; —dt, d > ma. 236 

In the opposite case we have to use a discrete sum because of a formal infrared divergence of the integral. Again, 
keeping only the main term in the expansion of the exponent in Eq. I|229(l . exp(& Yli=i ~~ 1 ~ T £i=i we 
find 

*m = {^^ /d /[r(- + l)] a [C(-)Y /d } m al m J, a<d< ma, (237) 
l c a a 

where a^J = Y^i>o(J^j=i^ 'i) m >d I ' £ T '■ (^ or ^ ne sa ^- e °^ simplicity, as in Eq. I|248|) . we again omit here a discussion of 
an obvious logarithmic factor that suppresses the ultraviolet divergence in the latter sum X)i>o ^ or ^ ne marginal case 
d = ma; see, e.g., Eq. I|250|) .) 

We conclude that all cumulants up to the order m < d/a have normal behavior, n m oc N , but for the higher orders, 
m > d/a, they acquire an anomalous growth in the thermodynamic limit, n m ~ k m oc ]\[ m < 7 / d . This result provides a 
simple classification of the relative strengths of the higher order fluctuation properties of the condensate in different 
traps. In particular, it makes it obvious that for all power-law traps with 1 < d/a < 2 the condensate fluctuations 
are not Gaussian, since 

oc N a -> const ^ at N -> oo, m > 3, (238) 

so that the asymmetry coefficient, 71 = ((no — n,o) 3 )/((no — "-o) 2 ) 3 ^ 2 7^ 0, and the excess coefficient, 72 = ((no — 
no) 4 )/((no — no) 2 ) 2 — 3 7^ 0, are not zero. Traps with d/a > 2 show Gaussian condensate fluctuations, since all 
higher-order cumulant coefficients K m /ft™^ 2 vanish, namely, 

oc N 1 -™/ 2 -> at N^oo if 3 < m < -, (239) 

a 



1/2 



m/2 
K 2 



N m ^-h)^Q a t N —>■ 00 if m>-. (240) 

(T 
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(For the sake of simplicity, we omit here an analysis of the special cases when d/a is an integer. It also can be 
done straightforwardly on the basis of the result l|229[) .') Very likely, a weak interaction also violates this non-robust 
property and makes properties of the condensate fluctuations in the traps with a relatively high dimension of the 
space, d/a > 2, similar to that of the box with the homogeneous Bose gas (see Section VI below), as it is stated below 
for the particular case of the harmonic traps. 

For traps with a space dimension lower than the critical value, d < er, it is known that a BEC phase transition 
does not exist (see, e.g., ^^). Nevertheless, even in this case there still exists a well-peaked probability distribution 
po(no) at low enough temperatures, so that the condition H21U|) is satisfied and our general result l|223|) describes this 
effect as well. In this case there is a formal infrared divergence in the corresponding integrals for all cumulants (|223[) . 
starting with the mean value. Hence, all of them should be calculated as discrete sums. For moderate temperatures 
we find approximately 

i) E(E^l) -(^-^EG^r) > ( 241 ) 

1>0 j=l 3=1 J 

so that higher cumulants have larger values. In particular, the mean number of non-condensed atoms is of the same 
order as the variance, 

KlS fi S N- no - lE(E^l)" 1 ~ E 

1>0 3=1 j=l J 



\ 



EG")' (^2) 



Therefore, until 



d _j d 

An2«iV, i.e., T<t:T c = N/J2(j2 hLU ^) ^ N /T,J— > ( 243 ) 

i>o j=i j=i J 



there is a well-peaked condensate distribution, y/An^ <C n . 

The marginal case of a trap with the critical space dimension, d = a, is also described by our result (EHJ), but we 
omit its discussion in the present paper. We mention only that there is also a formal infrared divergence, in this case 
a logarithmic divergence, and, at the same time, it is necessary to keep the whole exponent in Eq. ({229(1 : because 
otherwise in an approximation like i|241l) there appears an ultraviolet divergence. The physical result is that in such 
traps, e.g., in a one-dimensional harmonic trap, a quasi-condensation of the ideal Bose gas takes place at the critical 
temperature [UHl] T c ~ tvu)iN/ h\N . 

For very low temperatures, T -C £i, the second and higher energy levels in the trap are not thermally excited 
and atoms in the ideal Bose gas in any trap go to the ground level with an exponential accuracy. This situation is 
also described by Eq. I|229(l that yields fio ~ N and proves that all cumulants of the number-of-non-condensed-atom 
distribution become the same, since all higher-order generating cumulants exponentially vanish faster than Kx, 

d 

n m = £ m ~(m-l)!^ e - m ^'/ T , T« £l . (244) 

i=i 

The conclusion is that for the ideal Bose gas in any trap the distribution of the number of non-condensed atoms 
becomes Poissonian at very low temperatures. It means that the distribution of the number of condensed atoms is 
not Poissonian, but a "mirror" image of the Poisson's distribution. We see, again, that the complementary number of 
non-condensed atoms, n = N — hq, is more convenient for the characterization of the condensate statistics. A physical 
reason for this is that the non-condensed atoms in different excited levels fluctuate independently. 

The formulas l|229|) - (|244|l for the power-law traps contain all corresponding formulas for a box (d = 3, a = 2) 
and for a harmonic trap (d = 3, a = 1) in a 3-dimcnsional space as particular cases. It is worth to stress that 
BEC fluctuations in the ideal gas for the latter two cases are very different. In the box, if the temperature is larger 
than the trap energy gap, £\ <C T < T c , all cumulants, starting with the variance, m > 2, are anomalously large and 
dominated by the lowest energy modes, i.e., formally infrared divergent, 

Km « k m oc (77T c ) m 7V 2m / 3 , m > 2. (245) 
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Only the mean number of condensed atoms, 

/ / , ^/9\ 2ith 2 [ N \2/3 

n = N-n = N- Kl =N(l-(T/T c f/*), r c = _(__) , (246) 

can be calculated correctly via replacement of the discrete sum in Eq. (|223(l by an integral J °° ...d 3 k. The correct 
value of the squared variance, 

^<K-n ) 2 >=iV 4 / 3 (!) ^3^4/3 , S4 = g^ = 16.53, (247) 

can be calculated from Eq. 1|223[) only as a discrete sum. Thus, for the box the condensate fluctuations are anomalous 
and non-Gaussian even in the thermodynamic limit. To the contrary, for the harmonic trap with temperature much 
larger than the energy gap, E\ <C T < T c , the condensate fluctuations are Gaussian in the thermodynamic limit. 
This is because, contrary to the case of the homogeneous gas, in the harmonic trap only the third and higher order 
cumulants, m > 3, are lowest-energy-mode dominated, i.e., formally infrared divergent, 

1>0 



K ro «K TO «(m-l)!(!) m ^(wl) m cx7V™/ 3 , m> 3, (248) 



• T > 

1>0 

and they are small compared with an appropriate power of the variance squared 



C(3)Vr ( 



K2 =< (n - no) 2 >=^(^rfN. (249) 



The asymmetry coefficient behaves as 



_ k 3 _ <(n -n ) 3 > logiV 

and all higher normalized cumulants (m > 3) vanish in the thermodynamic limit, N — + oo, as follows: 

-TtL-O- (251) 



m/2 ]\[m/6 

It is important to realize that a weak interaction violates this non-robust property and makes properties of the 
condensate fluctuations in a harmonic trap similar to that of the homogeneous gas in the box (see Section VI below). 
For the variance, the last fact was first pointed out in |78| . 



D. Equivalent formulation in terms of the poles of the generalized Zeta function 

Cumulants of the BEC fluctuations in the ideal Bose gas, Eq. (|223|l . can be written in an equivalent form which 
is quite interesting mathematically (see j9|| and references therein). Namely, starting with the cumulant generating 
function lnH ea: (/3, z), where (3 = 1/kgT and z = e°^, 

1HS..G8, *) = - f>(l - *«p[-/J(e* ~ e )]) = f) £ e ^~ n ^ Z £ ^ , (252) 

v—\ y—1 n — 1 

we use the Mellin-Barnes transform 

er a = — / dta-*r(t) 

2™ J T - ioo 
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to write 



°° °° z n 1 r T +i°° J 

b3 - (ftz)= SS»^t to dtr(t) »^)i 



T+ioo 00 00 

\ " 1 \" 



1 

27r *ir^oo W ^ - £ 0)]* ^ 



*r(*)>.7S7i — —>,—■ (253) 



Recalling the series representation of the Bose functions g a (z) — J^^-i z n /»" and introducing the generalized, 
"spectral" Zeta function Z([3,t) — \fi{e v — eo)] _ *j we arrive at the convenient (and exact) integral representation 



1 



T + lOO 



lnS e:c (/3,z) = — / *r(t)ZC9,t) ffH .i(«). (254) 

Utilizing the well-known relations z-4-g a (z) = 3 Q -i( 2: ) and <? a (l) = C( a )i where C( a ) denotes the original Riemann 
Zeta function, Eq. (|254|) may be written in the appealing compact formula 



r+ioo 

dtT(t)Z(p,t)((t + l - k). (255) 



Thus, by means of the residue theorem, Eq. I|255|) links all cumulants of the canonical distribution in the condensate 
regime to the poles of the generalized Zeta function Z(/3,t), which embodies all the system's properties, and to the 
pole of a system-independent Riemann Zeta function, the location of which depends on the order k of the respective 
cumulant. The formula (|255|1 provides a systematic asymptotic expansion of the cumulants /tfe(/3) through the residues 
of the analytically continued integrands, taken from right to left. The large-system behavior is extracted from the 
leading pole, finite-size corrections are encoded in the next-to-leading poles, and the non-Gaussian nature of the 
condensate fluctuations is definitely seen. The details and examples of such analysis can be found in [9fil |. 

Concluding the Section V, it is worthwhile to mention that previously only first two moments, K\ and «2, were 
analyzed for the ideal gasjUlilililllllllllillill^, and the known results coincide with ours. Our explicit 
formulas provide a complete answer to the problem of all higher moments of the condensate fluctuations in the ideal 
gas. The canonical-ensemble quasiparticle approach, taken together with the master equation approach gives a fairly 
complete picture of the central moments. 

For the interacting Bose gas this problem becomes much more involved. We address it in the next, last part of this 
review within a simple approximation that takes into account one of the main effects of the interaction, namely, the 
Bogoliubov coupling. 



VI. WHY CONDENSATE FLUCTUATIONS IN THE INTERACTING BOSE GAS ARE 
ANOMALOUSLY LARGE, NON-GAUSSIAN, AND GOVERNED BY UNIVERSAL INFRARED 

SINGULARITIES ? 

In this Section, following the Refs. pol. l2lj. the analytical formulas for the statistics, in particular, for the char- 
acteristic function and all cumulants, of the Bose-Einstein condensate in the dilute, weakly interacting gases in the 
canonical ensemble are derived using the canonical ensemble quasiparticle method. We prove that the ground-state 
occupation statistics is not Gaussian even in the thermodynamic limit. We calculate the effect of Bogoliubov cou- 
pling on suppression of ground-state occupation fluctuations and show that they are governed by a pair-correlation, 
squeezing mechanism. 

It is shown that the result of Giorgini, Pitaevskii, and Stringari (GPS) [78j for the variance of condensate fluctuations 
is correct, and the criticism of Idziaszek and others 97, 98] is incorrect. A crossover between the interacting and ideal 
Bose gases is described. In particular, it is demonstrated that the squared variance of the condensate fluctuations 
for the interacting Bose gas, Eq. (|271|l . tends to a half of that for the ideal Bose gas, Eq. I|247|) . because the atoms 
are coupled in strongly correlated pairs such that the number of independent degrees of freedom contributing to the 
fluctuations of the total number of excited atoms is only 1/2 the atom number N . This pair correlation mechanism 
is a consequence of two-mode squeezing due to Bogoliubov coupling between k and — k modes. Hence, the fact that 
the fluctuation in the interacting Bose gas is 1/2 of that in the ideal Bose gas is not an accident, contrary to the 
conclusion of GPS. Thus, there is a deep (not accidental) parallel between the fluctuations of ideal and interacting 
bosons. 
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Finally, physics and universality of the anomalies and infrared singularities of the order parameter fluctuations for 
different systems with a long range order below a critical temperature of a second order phase transition, including 
strongly interacting superfluids and ferromagnets, is discussed. In particular, an effective nonlinear a model for the 
systems with a broken continuous symmetry is outlined and the crucial role of the Goldstone modes fluctuations 
combined with an inevitable geometrical coupling between longitudinal and transverse order parameter fluctuations 
and susceptibilities in the constrained systems is demonstrated. 



A. Canonical-ensemble quasiparticles in the atom-number-conserving Bogoliubov approximation 



We consider a dilute homogeneous Bose gas with a weak interatomic scattering described by the well-known Hamil- 
tonian 

h 2 k 2 1 

H = X, ~2M h ^ + W ^ ( k 3k4|^|kik 2 )a+ a+ a k2 a kl , (256) 

k {k ; } 

where V = I? is a volume of a box confining the gas with periodic boundary conditions. The main effect of the 
weak interaction is the Bogoliubov coupling between bare canonical-ensemble quasiparticles, /3k = /3^~ak, via the 
condensate. It may be described, to a first approximation, by a quadratic part of the Hamiltonian (|256Jl . i.e., by the 
atom-number-conserving Bogoliubov Hamiltonian 93] 

Hb _ + + (ft. + 'WW. )^ + ^fc ^iTgHjfe + „,.), ( 257 ) 

where we will make an approximation % ~ no > 1, which is consistant with our main assumption (|210(l of the 
existence of a well peaked condensate distribution function. Then, the Bogoliubov canonical transformation, 

1 Au 
/?k = "kfrk + Vkb\; «k = —i===, v k 



describes the condensate canonical-ensemble quasiparticles which have a "gapless" Bogoliubov energy spectrum and 
fluctuate independently in the approximation l|257|) . since 

k^O 

In other words, we again have an ideal Bose gas although now it consists of the dressed quasiparticles which are 
different both from the atoms and bare (canonical-ensemble) quasiparticles introduced in Section V. Hence, the 
analysis of fluctuations can be carried out in a similar fashion to the case of the non-interacting, ideal Bose gas of 
atoms. This results in a physically transparent and analytical theory of BEC fluctuations that was suggested and 
developed in ^H^- 

The only difference with the ideal gas is that now the equilibrium density matrix, 

P = n M o/3 k! Pk - e~ s ^/ T (l - e^/ T ) , (260) 

is not diagonal in the bare atomic occupation numbers, the statistics of which we are going to calculate. This feature 
results in the well-known quantum optics effect of squeezing of the fluctuations. The number of atoms with coupled 
momenta k and — k is determined by the Bogoliubov coupling coefficients according to the following equation: 

ajak + a± k a- k = /3^/3 k + /3+ k /3_ k = (u£ + vl)(b+b k + b±J>-k) + 2u k v k (b+b\ + bJ)-k) + 2v£. (261) 
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B. Characteristic function and all cumulants of BEC fluctuations 



The characteristic function for the total number of atoms in the two, k and k, modes squeezed by Bogoliubov 
coupling is calculated in [20I l2l| as 

e ±k (u) = Tr(e %u{ ^ $k+ ^- k) e- £k{i i' bk+i -^- k)/T (l ~ e- £k/T 



(z(A k )-l)(*(-A k )-l) A k -e°*/ T 

< A *) = T „,ut ; ■ ( 262 ) 



(z(A k )-e™){z(-A k )-e™)> K A k e^/ T - 1' 

The term "squeezing" originates from the studies of a noise reduction in quantum optics (see the discussion after 

E q . arm ). 

The characteristic function for the distribution of the total number of the excited atoms is equal to the product of 
the coupled-mode characteristic functions, O n (u) = n k ^o.mo!i{±k}0±k('«), since different pairs of (k, — k)-modes 
are independent to the first approximation i|257[l . The product II runs over all different pairs of (k, — k)-modes. 

It is worth noting that by doing all calculations via the canonical-ensemble quasiparticles (Section V) we automat- 
ically take into account all correlations introduced by the canonical-ensemble constraint. As a result, similar to the 
ideal gas (Eq. (|223|l % ). we obtain the explicit formula for all cumulants in the dilute weakly interacting Bose gas, 

1 11 r 

km = 2 (m_1)! 5W)-ir + w-^-iH' Kr = E 4m) ~ Km - (263) 

In comparison with the ideal Bose gas, Eq. I|223[) . we have effectively a mixture of two species of atom pairs with 
z(±A k ) instead of exp(£k/T). 

It is important to emphasize that the first equation in l|263|) . m = 1, is a non-linear self-consistency equation, 



N - n = Kl (n ) ^ £ _ * - (264) 

to be solved for the mean number of ground-state atoms fio(T), since the Bogoliubov coupling coefficient (|258J) . A k , 
and the energy spectrum (|259|l . £k, are themselves functions of the mean value n®. Then, all the other equations in 
(I263|l . m > 2, are nothing else but explicit expressions for all cumulants, m > 2, if one substitutes the solution of 
the self-consistency equation <|264[) for the mean value n . The Eq. (|264|) . obtained here for the interacting Bose gas 
(|257|l in the canonical-ensemble quasiparticle approach, coincides precisely with the self-consistency equation for the 
grand-canonical dilute gas in the so-called first-order Popov approximation (see a review in |36j). The latter is well 
established as a reasonable first approximation for the analysis of the finite-temperature properties of the dilute Bose 
gas and is not valid only in a very small interval near T c , given by T c — T < a(N /V) l '^T c -C T c , where a — MUo/Anh 2 
is a usual s-wave scattering length. The analysis of the Eq. (|264|) shows that in the dilute gas the self-consistent value 
no(T') is close to that given by the ideal gas model, Eq l|246l) . and for very low temperatures goes smoothly to the 
value given by the standard Bogoliubov theory j^O, l94l| for a small condensate depletion, N — no <C N. This is 
illustrated by Fig. 114b . (Of course, near the critical temperature T c the number of excited quasiparticles is relatively 
large, so that along with the Bogoliubov coupling 1)2571) other, higher order effects of interaction should be taken into 
account to get a complete theory.) Note that the effect of a weak interaction on the condensate fluctuations is very 
significant (see Fig. I14b -d) even if the mean number of condensed atoms changes by relatively small amount. 



C. Surprises: BEC fluctuations are anomalously large and non-Gaussian even in the thermodynamic limit 

According to the standard textbooks on statistical physics, e.g., [3(], 0, H3' any extensive variable C of a 
thermodynamic system has vanishing relative root-mean-square fluctuations. Namely, in the thermodynamic limit, a 
relative squared variance ((C — C) 2 ) / C 2 — (C 2 ) / (C) 2 — 1 oc V~ 1 goes inversely proportional to the system volume V, 
or total number of particles N. This fundamental property originates from the presence of a finite correlation length 
£ that allows us to partition a large system into an extensive number V/£ 3 of statistically independent subvolumes, 
with a finite variance in each subvolume. As a result, the central limit theorem of probability theory yields a Gaussian 

distribution for the variable C with a standard scaling for variance, y {{C — C) 2 ) ~ V 1 ^ 2 . Possible deviations from 
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FIG. 14: Temperature scaling ol the first four cumulants, the mean value no/N — N — ni/N, the variance 
{(n - n ) 2 ) 1/2 /N 1/2 , the third central moment -k}/ 3 /N 1/2 = ((no - n ) 3 ) 1/3 /N 1/2 , the fourth cumulant \k 4 \ 1/4 /N 1/2 = 
| ((no — no) 4 ) — 3k 2 \/N 2 , of the ground-state occupation fluctuations for the dilute weakly interacting Bose gas (Eq. 1)263^ . 
with UoN 1 / 3 1 'e\V — 0.05 (thick solid lines), as compared with Eq. 1223H (thin solid lines) and with the exact recursion relation 
1801 (dot-dashed lines) for the ideal gas in the box; N = 1000. For the ideal gas our results (thin solid lines) are almost 
indistinguishable from the exact recursion calculations (dot-dashed lines) in the condensed region, T < T C (N). Temperature is 
normalized by the standard thermodynamic-limit critical value T C (N = oo) that differs from the finite-size value T C (N), as is 
clearly seen in graphs. 
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this general rule are especially interesting. It turns out that BEC in a trap is one of the examples of such peculiar 
systems. A physical reason for the anomalously large and non-Gaussian BEC fluctuations is the existence of the long 
range order below the critical temperature of the second order phase transition. 

Let us show in detail that the result (|263[1 implies, similar to the case of the ideal homogeneous gas (Section V), 
that the ground-state occupation fluctuations in the weakly interacting gas are not Gaussian in the thermodynamic 
limit, and anomalously large. The main fact is that the anomalous contribution to the BEC fluctuation cumulants 
comes from the modes which have the most negative Bogoliubov coupling coefficient « — 1 since in this case one 
has z{A]s) — > 1, so that this function produces a singularity in the first term in Eq. (|263[1 . (The second term in 
Eq. (|263(l never makes a singular contribution since always z{— A-^) < —I or z{— A^) > exp(£k/"T).) These modes, 
with Ak ss —I, exist only if the interaction energy g — uqU^/V in the Bogoliubov Hamiltonian (|257[) is larger than 

the energy gap between the ground state and the first excited states in the trap, e\ — 
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(265) 



and are the infrared modes with the longest wavelength of the trap energy spectrum. In terms of the scattering 
length for atom-atom collisions, a = MUk/Airh 2 , the latter condition l|265[l coincides with a familiar condition for the 
Thomas-Fermi regime, in which the interaction energy is much larger then the atom's kinetic energy. 
Let us use a representation which is obvious from Eqs. I|258|l and l|259|) . 



z(A k ) - 1 



2ei, 



.9(1 - ^k) 



A(l - A v ) 
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— tanhf — ^ 
9 V2T 



(266) 



where in the last, approximate equality we set A^ ~ — 1, and neglect the contribution from the second term in Eq. 
(|263|l . assuming that the singular contribution from the modes with A^ w — 1 via the first term in Eq. I|263[l is 
dominant. Then, for all infrared-dominated cumulants of higher orders to > 2, the result (|263|l reduces to a very 
transparent form 
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(267) 



Finally, using the Bogoliubov energy spectrum in Eq. (|259|l . e k = E\ ^/l 4 + (2g/ei)l 2 with a set of integers 1 = 
{Ixi ly,lz), we arrive to the following simple formulas for the higher order generating cumulants in the Thomas-Fermi 
regime (126511 : 
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(269) 
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(270) 



for high, moderate, and very low temperatures T respectively, as compared to the geometrical mean of the interaction 
and gap energies in the trap, yj g£i/2. 

In particular, for the Thomas-Fermi regime (|265J) and relatively high temperatures, the squared variance, as given 
I'.v K(|. 
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(271) 



scales as ((no — no) 2 ) cx iV 4 / 3 . Here the arrow indicates the limit of sufficiently strong interaction, g 3> E\ and 

T># > E X . 
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The behavior l|271|) , l|245[) , and (|275() is essentially different from that of the normal fluctuations of most extensive 
physical observables, which are Gaussian with the squared variance proportional to N. The only exception for the 
Eqs. (j267|l - l|27Q(l is the low temperature limit of the variance that is not infrared-dominated and should be calculated 
not via the Eq. (|27U[) . but directly from Eq. (|263[) using the fact that all modes are very poorly occupied at low 
temperatures, i.e., exp(— s-^/T) -C 1, if g ^> 2T 2 /ei,ei. Thus, we immediately find 

The above results extend and confirm the result of the pioneering paper w here only the second m omen t, 
((n — n ) 2 ), was calculated. (The result of [71] was rederived by a different way in |100| . and generalized in |lOl| .1 
The higher-order cumulants K m , m > 2, given by Eqs. I|263[) and (|267(l - l|270|l . are not zero, do not go to zero in the 
thermodynamic limit and, moreover, are relatively large compared with the corresponding exponent of the variance 
( Km ) m / 2 that proves and measures the non-Gaussian character of the condensate fluctuations. For the Thomas-Fermi 
regime (|265|) and relatively high temperatures, the relative values of the higher order cumulants are given by Eq. 
(USB} as 

-1,!^, ("3) 

where 1 = (l x ,l yi l z ) are integers, S4 ~ 16.53, sq w 8.40, and S2m ~ 6 for m 1. In particular, the asymmetry 
coefficient of the ground-state occupation probability distribution 

71 = ((n - n Q ) 3 )/((n - n ) 2 f' 2 = -k 3 /k 3 2 /2 » 2V2s 6 /( S4 ) 3/2 « 0.35 (274) 
is not very small at all. 

This non-Gaussian statistics stems from an infrared singularity that exists for the fluctuation cumulants n m , m > 2, 
for the weakly interacting gas in a box despite of the acoustic (i.e., linear, like for the ideal gas in the harmonic 
trap) behavior of the Bogoliubov-Popov energy spectrum (|259|l in the infrared limit. The reason is that the excited 
mode squeezing (i.e., linear mixing of atomic modes) via Bogoliubov coupling affects the BEC statistics i|263|l for 
the interacting gas also directly (not only via a modification of the quasiparticle energy spectrum), namely via the 
function 12621) . z(Ay : ), which is different from a simple exponent of the bare energy, exp(ek/T), that enters into the 
corresponding formula for the non-interacting gas 12231) . 



D. Crossover between ideal and interaction-dominated BEC: quasiparticles squeezing and pair correlation 



Now we can use the analytical formula l|263[) to describe explicitly a crossover between the ideal-gas and interaction- 
dominated regimes of the BEC fluctuations. Obviously, if the interaction energy g = fioUk/V is less than the energy 
gap between the ground state and the first excited state in the trap, g < e\, the Bogoliubov coupling (|258|l becomes 
small for all modes, |^4k| "C 1, so that both terms in Eq. (|263|l give similar contributions and all fluctuation cumulants 
k m tend to their ideal gas values in the limit of vanishing interaction, g <C £1. Namely, in the near-ideal gas regime 
fioa/L -C 1 the squared variance linearly decreases from its ideal-gas value with an increase of the weak interaction 
as follows: 
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With a further increase of the interaction energy over the energy gap in the trap, g > e±, the essential differences 
between the weakly interacting and ideal gases appear. First, the energy gap is increased by the interaction, that is 



£1 



yjel + 2e 1 n (T)Uo/V > e x = (^) * JL , 



(276) 



so that the border T ~ £1 between the moderate temperature and very low temperature regimes is shifted to a 
higher temperature, T ~ y/ gei/2. Thus, the interaction strength g determines also the temperature T ~ y ge\j2 
above which another important effect of the weak interaction comes into play. Namely, according to Eqs. I|267|) - (|270|l . 
the suppression of all condensate- fluctuation cumulants by a factor of 1/2, compared with the ideal gas values (see 
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FIG. 15: The probability distribution l|278|l P2 as a lunction of the number of atoms in k and — k modes, rik + n_k, for the 
interaction energy Uono/V = 10 3 £k and temperature T — e^. The pair correlation effect due to Bogoliubov coupling in the 
weakly interacting Bose gas is clearly seen for low occupation numbers, i.e., even occupation numbers nk + n_k are more 
probable than odd ones. 



Eq. I|247|) ). takes place for moderate temperatures, £\ <C T < T C) when a strong coupling (A k 
dominates in Eq. (|263|l . The factor 1/2 comes from the fact that, according to Eq. (|262|l . 
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so that the first term in Eq. I|263(l is resonantly large but the second term is relatively small. In this case, the effective 
energy spectrum, which can be introduced for the purpose of comparison with the ideal gas formula l|223l) , is 



= T]n(z(A k ))~±(l + A*)e k 



\elv/u Q n * A_ 



(277) 



That is, the occupation of a pair of strongly coupled modes in the weakly interacting gas can be characterized by 
the same effective energy spectrum as that of a free atom. It is necessary to emphasize that the effective energy 
e^ff = Tln(z(^4 k )), introduced in Eq. 1)277(1 . describes only the occupation of a pair of bare atom excitations k and 
— k (see Eqs. <|260[) - (|262ll > ) and, according to Eq. I|263|) . the ground-state occupation. That is, it would be wrong to 
reduce the analysis of the thermodynamics and, in particular, the entropy of the interacting gas to this effective energy. 
Thermodynamics is determined by the original energy spectrum of the dressed canonical-ensemble quasiparticles, Eq. 



This remarkable property explains why the ground-state occupation fluctuations in the interacting gas in this 
case are anomalously large to the same extent as in the non-interacting gas except factor of 1/2 suppression in the 
cumulants of all orders. These facts were considered in to be an accidental coincidence. We see now that, roughly 
speaking, this is so because the atoms are coupled in strongly correlated pairs such that the number of independent 
stochastic occupation variables ( "degrees of freedom" ) contributing to the fluctuations of the total number of excited 
atoms is only 1/2 the atom number AT. This strong pair correlation effect is clearly seen in the probability distribution 
of the total number of atoms in the two coupled k and — k modes, 
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(see Fig. 1150 . The latter formula follows from Eq. I|262|) . Obviously, a higher probability for even occupation numbers 
ri k + n_ k as compared to odd numbers at low occupations means that the atoms in the k and — k modes have a 
tendency to appear or disappear simultaneously, i.e., in pairs. This is a particular feature of the well-studied in 
quantum optics phenomenon of two-mode squeezing (see, e.g., |l8l \3l\). This squeezing means a reduction in the 
fluctuations of the population difference ri k — H- k , and of the relative phases or so-called quadrature-phase amplitudes 
of an interacting state of two bare modes a k and a_ k compared with their appropriate uncoupled state, e.g., coherent 
or vacuum state. The squeezing is due to the quantum correlations which build up in the bare excited modes via 
Bogoliubov coupling l|258|l and is very similar to the noise squeezing in a non-degenerate parametric amplifier studied 
in great details by many authors in quantum optics in 80s 0, y^] ■ We note that fluctuations of individual bare 
excited modes are not squeezed, but there is a high degree of correlation between occupation numbers in each mode. 
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It is very likely that in the general case of an arbitrary power-law trap the interaction also results in anomalously large 
fluctuations of the number of ground-state atoms, and a formal infrared divergence due to excited mode squeezing via 
Bogoliubov coupling and renormalization of the energy spectrum. In the particular case of the isotropic harmonic trap 
this was demonstrated in [T^ for the variance of the condensate fluctuations. Therefore, the ideal gas model for traps 
with a low spectral index a < d/2 (such as a three-dimensional harmonic trap where a = 1 < d/2 = 3/2), showing 
Gaussian, normal thermodynamic condensate fluctuations with the squared variance proportional to N instead of 
anomalously large fluctuations (see Eqs. I|232|) . 1234|) and l|239|) . I|240|) 'l. is not robust with respect to the introduction 
of a weak interatomic interaction. 

At the same time, the ideal gas model for traps with a high spectral index a > d/2 (e.g., for a three-dimensional 
box with a = 2 > d/2 = 3/2) exhibits non-Gaussian, anomalously large ground-state occupation fluctuations with a 
squared variance proportional to N 2<I / d N (see Eqs. (|233|l . (|234|l % ) similar to those found for the interacting gas. 
Fluctuations in the ideal Bose gas and in the Bogoliubov Bose gas differ by a factor of the order of 1, which, of course, 
depends on the trap potential and is equal to 1/2 in the particular case of the box, where An,Q oc N 4 ' 3 . We conclude 
that, contrary to the interpretation formulated in [x^, similar behavior of the condensate fluctuations in the ideal 
and interacting Bose gases in a box is not accidental, but is a general rule for all traps with a high spectral index 
(7 > d/2, or a relatively low dimension of space, d < 2a. 

As follows from Eq. (|263|l . the interaction essentially modifies the condensate fluctuations also at very low tem- 
peratures, T <C £i (see Figs. H~4Tl . Namely, in the interacting Bose gas a temperature-independent quantum noise, 



K m -> R m (T = 0)^0, m > 2, (279) 

additional to vanishing (at T — > 0) in the ideal Bose gas noise, appears due to quantum fluctuations of the excited 
atoms, which are forced by the interaction to occupy the excited levels even at T = 0, so that n^(T = 0) ^ 0. Thus, 
in the limit of very low temperatures the results of the ideal gas model (Section V) are essentially modified by weak 
interaction and do not describe condensate statistics in the realistic weakly interacting Bose gases. 

The temperature scaling of the condensate fluctuations described above is depicted in Fig. [21 both for the weakly 
interacting and ideal gases. A comparison with the corresponding quantities calculated numerically from the exact 
recursion relation in Eqs. I|79|l and (|80|l for the ideal gas in a box is also indicated. It is in good agreement with 
our approximate analytical formula (|263|l for all temperatures in the condensed phase, T < T c , except of a region 
near to the critical temperature, T w T c . It is worth stressing that the large deviations of the asymmetry coefficient, 
7i = ((no — n o) 3 )/((no — n o) 2 ) 3 ^ 2 , and of the excess coefficient, 72 = ((no — n ) 4 } / ((n — no) 2 ) 2 — 3 from zero, which 
are of the order of 1 at T ~ T c /2 or even more at T as and T « T C) indicate how far the ground-state occupation 
fluctuations are from being Gaussian. (In the theory of turbulence, the coefficients 71 and 72 are named as skewness 
and flatness, respectively.) This essentially non-Gaussian behavior of the ground-state occupation fluctuations remains 
even in the thermodynamic limit. 

Mesoscopic effects near the critical temperature are also clearly seen in Fig. ^[ and for the ideal gas are taken into 
account exactly by the recursion relations 17911 , (181) I) . The analytical formulas (|263[) take them into account only via a 
finite-size effect, L oc TV -1 / 3 , of the discreteness of the single-particle spectrum e^. This finite-size (discreteness) effect 
produces, in particular, some shift of the characteristic BEC critical temperature compared to its thermodynamic- 
limit value, T c . For the case shown in Fig. 1141 it is increased by a few per cent. Similarly to the mean-number "grand" 
canonical approximation described at the end of Section III.B for the case of the ideal gas, the canonical-ensemble 
quasiparticle approach can partially accommodates for the mesoscopic effects by means of the grand-canonical shift of 
all quasiparticle energies by a chemical potential, £k = £k — M- I n this case the self-consistency equation (|264[) acquires 
an additional nonlinear contribution due to the relation exp(— 0fi) = 1 + 1/no- Obviously, this "grand" canonical 
approximation works only for T < T c , whereas at T > T c we can use the standard grand canonical approach, since the 
ground-state occupation is not macroscopic above the critical temperature. However, the "grand" canonical approach 
takes care only of the mean number of condensed atoms no, and does not improve the results of the canonical-ensemble 
quasiparticle approach for BEC fluctuations. 



E. Universal anomalies and infrared singularities of the order parameter fluctuations in the systems with a 

broken continuous symmetry 



The result that there are singularities in the central moments of the condensate fluctuations, emphasized in 20, 21] 
and discussed above in detail for the BEC in a trap, can be generalized for other long-range ordered systems below 
the critical temperature of a second-order p hase trans ition, including strongly interacting systems. That universality 
of the infrared singularities was discussed in |lQCt Il0l| and can be traced bac k to a well known property of an infrared 
singularity in a longitudinal susceptibility x\\(k) °f such systems [3(], Il02j . The physics of these phenomena 
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is essentially determined by long wavelength phase fluctuations, which describe the non-condensate statistics, and 
is intimately related to the fluctuation-dissipation theorem, Bogoliubov's 1/fc 2 theorem for the static susceptibility 
of superfluids, and the presence of Goldstone modes, as will be detailed in the following. 



Long wavelength phase fluctuations, fluctuation-dissipation theorem, and Bogoliubov's 1/fc 2 theorem 



First, let us refer to a microscopic derivation given in [lj03j that demonstrates the fact that the phase fluctuations 
alone dominate the low-energy physics. Also, let us assume that, as was shown by Feynman, the spectrum of 
excitations in the infrared limit k — > is exhausted by phonon-like modes with a linear dispersion u>k = ck, where c 
is an actual velocity of sound. Then, following a textbook [30j, we can approximate an atomic field operator via a 
phase fluctuation operator as follows: 

tt(x) = v^e^W, 0(x) = K/2^ to( ) 1/2 ^(ffi)- 1/2 (c k e lkx + c+er^), (280) 

k^O 

where h is the bare condensate density, n tot the mean particle density, m the particle mass, Ck and are phonon 
annihilation and creation operators, respectively. An omission of the k — term in the sum in Eq. I|28U|) can be 
rigorously justified on the basis of the canonical-ensemble quasiparticle approach (see Section V) and is related to the 
fact that a global phase factor is irrelevant to any observable gauge-invariant quantity. In a homogeneous supcrfluid, 
the renormalized condensate density is determined by the long range order parameter, 

no = Hm <* + (Ax)*(0)) w h Q exp(-(^ 2 (0))), (281) 

Ax — >oo 

if we use the approximation (|28()f> and neglect phonon interaction. Thus, the mean number of condensate particles 
is depicted with an increase of temperature in accordance with an increase of the variance of the phase fluctuations, 
no(T) oc exp(— (<^ 2 (0))t), which yields a standard formula |30( for the thermal depletion, fio(T) — no(T = 0) = 
-n (T = Q)m(k B T) 2 /12n tot ch 3 . 

Similar to the mean value, higher moments of the condensate fluctuations of the homogeneous superfluid in the 
canonical ensemble can be represented in terms of the phase fluctuations of the non-condensate via the operator 

h = N/V - n = J Ai>+{x)Ai>(x)d 3 x, (282) 

where 

Mf{x) = 9(x) - V^Tj = v^(e**(aO - e -^ 2 (°)>/ 2 ). (283) 

In particular, the vari ance of the condensate occupation is determined by the correlation function of the phase 
fluctuations as follows |100| : 

((no - rlo) 2 ) w 2n 2 e- 2 ^ 2 W> J {<p{x)<f{y)) 2 Sx<fy. (284) 
The last formula can be evaluated with the help of a classical form of the fluctuation-dissipation theorem, 

(0(x)0(y)) 2 d 3 xd 3 y = ^(fc s T Xw (k)) 2 . (285) 

For a homogeneous superfluid, the static susceptibility in the infrared limit k — > does not depend on the intera ction 
strength and, in accordance with Bogoliubov's l/fc 2 -theorem for the static susceptibility of superfluids, is equal to |l00j 



Xw (k) = m/n tot h 2 k 2 , k -> 0. (286) 

Then Eq. I|284|) immediately yields the squared variance of the condensate fluctuations in the superfluid with arbitrary 
strong interaction in exactly the same form as is indicated by the arrow in Eq. (|271|l . only with the additional factor 
(no(0) / ntot) 2 ■ The latter factor is almost 1 for the BEC in dilute Bose gases but can be much less in liquids, for 
example in A He superfluid it is about 0.1. Thus, indeed, the condensate fluctuations at low temperatures can be 
calculated via the long wavelength phase fluctuations of the non-condensate. 
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2. Effective nonlinear a model, Goldstone modes, and universality of the infrared anomalies 



Following |lOl| . let us use an effective nonlinear a model |l04j to demonstrate that the infrared singularities and 
anomalies of the order parameter fluctuations exist in all systems with a broken continuous symmetry, independently 
on the interaction strength, and are similar to that of the BEC fluctuations in the Bose gas. That model describes 
directional fluctuations of an order parameter ^(x) 



,(°) 



Q(x) with a fixed magnitude mi°' in terms of an Nq 



component unit vector Q(x). It is inspired by the classical theory of spontaneous magnetization in ferromagnets. The 
constraint = 1 suggests a standard decomposition of the order parameter 



n(ar) = {no(x),SU(x);i = l,...,N a - 1}, Cl (x) 



\ 



JVn-l 

i=i 



(287) 



into a longitudinal component ilo(x) and Nn — 1 transverse Goldstone fields Cli(x). The above a model constraint, 
£Iq(x) = 1 — J2i resembles the particle- number constraint no = N — X)k^o " k m -^q- tUP f° r the many-body 

atomic Bose gas in a trap. Although the former is a local, more stringent, constraint and the latter is only a global, 
integral constraint, in both cases it results in the infrared singularities, anomalies, and non-Gaussian properties of the 
order parameter fluctuations. Obviously, in the particular case of a homogeneous system the difference between the 
local and integral constraints disappears entirely. Within the a model, a superfluid can be described as a particular 
case of an Nq — 2 system, with the superfluid (condensate) and normal (non-condensate) component. 

At zero external field, the effective action for the fluctuations of the order parameter is 5[f2] = 
(p s /2T) J[\7il(x)] 2 d 3 x, where the spin stiffness p s is the only parameter. Below the critical temperature T c , a contin- 
uous symmetry becomes broken and there appears an intensive nonzero order parameter with a mean value 



■ $(0)}d 3 x = Vm\ -> Vm\ 



(288) 



which gives the spontaneous magnetization m 2 of the infinite system in the thermodynamic limit V — > oo. The 
leading long distance behavior of the two-point correlation function G(x) = (^(x) ■ ^(0)) may be obtained from a 
simple Gaussian spin wave calculation. Assuming low enough temperatures, we can neglect the spin wave interactions 
and consider the transverse Goldstone fields O, (k) as the Gaussian random functions of the momentum k with the 
correlation function (f2^(k)fV(k / )) = 6i,i>Su,-u'T/ p s k 2 . As a result, the zero external field correlation function below 
the critical temperature is split into longitudinal and transverse parts G(x) = m 2 s [\ + Gii (x) + (Afo — l)G±(x)], where 
m 2 = G(oo) now is the renormalized value of the spontaneous magnetization. To the lowest nontrivial order in the 
small fluctuations of the Goldstone fields, the transverse correlation function decays very slowly with a distance r, 
G± oc T/p s r, in accordance with Bogoliubov's l/fc 2 -theorem of the divergence of the transverse susceptibility in the 
infrared limit, 



X±(k) =m 2 s G ± (k)/T = m 2 /p s k 2 . 



(289) 



(290) 



The longitudinal correlation function is simply related to the transverse one |1£ 

G {] (x) « j£>2 (x) ^2 (0))c = l {Nn i )G l {xl 

and, hence, decays slowly with a 1/r 2 power law. That means that contrary to the naive mean field picture, where 
the longitudinal susceptibility X\\{k — * 0) below the critical temperat ure i s finite, the Eq. (|290|l leads to an infrared 
singularity in the longitudinal susceptibility and correlation function |102| . 



X\\(k^0)~T/p 2 k; 



Gii 



1/r 2 . 



(291) 



Although the Eq. (|290|l is obtained by means of perturbation theory, the result for t he sl ow power law decay of the 
longitudinal correlation function, Eq. Q291[) . holds for arbitrary temperatures T < T c |l04j . 
Knowing susceptibilities, it is immediately possible to find the variance of the operator 



M s =V~ 



$(x) ■ y(y)d 3 xd 3 y 



(292) 



V JV 



that describes the fluctuations of the spontaneous magnetization in a finite system at zero external field. Its mean 
value is given by Eq. (|288|l as (M s ) — M s = Vm 2 L . Its fluctuations are determined by the connected four-point 
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correlation function G4 = (^>(xi) ■ $ (2:2) ^(x^) ■ ^(x4)) c . In the 4-th order of the perturbation theory for an infinite 
system, it may be expressed via the squared transverse susceptibilities as follows: 

G 4 = \m 4 s (N n - - a*) - Gj_( Xl - x 4 ) - G±(x 2 - x 3 ) + G ± (x 2 - x 4 )] 2 . (293) 

For a finite system, it is more convenient to do similar calculations in the momentum representation and replace all 
integrals by the corresponding discrete sums, which yields the squared variance 

((M s - M s f) = 2m 4 (N n - l)T 2 p- 2 £ fc~ 4 oc TV 4 / 3 . (294) 

This expression has the same anomalously large scaling oc V 4 ! 3 in the thermodynamic limit V — > 00 (due to the same 
infrared singularity) as the variance of the BEC fluctuations in the ideal or interacting Bose gas in Eqs. I|247|) or (|271|l . 
respectively. Again, although the Eq. (|294|l is obtained by means of perturbation theory, the latter scaling is universal 
below the critical temperature, just like the l/k infrared singularity of the longitudinal susceptibility. Of course, very 
close to T c there is a crossover to the critical singularities, as discussed, e.g., in |l05]. Thus, the average fluctuation 

of the order parameter is still vanishing in the thermodynamic limit J ((M s — M s ) 2 ) / M s oc V~ x l 3 — > 0; that is, the 
order parameter (e.g., spontaneous magnetization or macroscopic wave function) is still a well defined self-averaging 
quantity. However, this self-averaging in systems with a broken continuous symmetry is much weaker than expected 
naively from the standard Einstein theory of the Gaussian fluctuations in macroscopic thermodynamics. Note that, 
in accordance with the well-known Hohenberg-Mermin- Wagner theorem, in systems with lower dimensions, d < 2, 
the strong fluctuations of the direction of the magnetization completely destroy the long range order, and the self- 
averaging order parameter does not exist anymore. An important point also is that the scaling result in Eq. I|294l) holds 
for any dynamics and temperature dependence of the average order parameter M S (T) although the function M S (T) 
is, of course, different, say, for ferromagnets, anti-ferromagnets, or a BEC in different traps. It is the constraint, either 
I SI I = 1 in the a model or N = ho + J^k^o " k m ^ ne BEC, that predetermines the anomalous scaling in Eq. (|294|) . The 
temperature dependence of the variance in Eq. I|294|) at low temperatures is also universal, since p s — > const at T — > 0. 
The reason for this fact is that the dominant finite size dependence is determined by the leading low energy constant 
in the effective field theory for fluctuations of the order parameter, which is precisely p s in the effective action S [Q]. 

3. Universal scaling of condensate fluctuations in superfluids 

In homogeneous superfluids the translational invariance requires the superfluid density to be equal to the full density 
n to t, so that the associated stiffness p s (T — > 0) = h 2 n to t/fn is independent of the interaction strength. Thus, Eq. I|294|) 
in accord with the Eq. (|271|l yields the remarkable conclusion that the relative variance of the ground-state occupation 
at low temperatures is a universal function of the density and the thermal wavelength At = /i/v / 27rmT, as well as the 
system size L = V 1 ! 3 and the boundary conditions, 

V((no-h )2)/n Q = B/(n tot \ 2 T L) 2 . (295) 

The boundary conditions determine the low-energy spectrum of quasiparticles in the trap in the infrared limit and, 
hence, the numerical prefactor B in the infrared singularity of the variance, namely, the coefficient B in the sum 
£k#o fc ~ 4 = BV 4/3 /8n 2 . In particular, in accord with Eqs. (|2T7|l and lf27Tl) . one has B = 0.8375 for the box with 
periodic boundary conditions, and B — %E 3 (2)/tt 2 ~ 0.501 for the box with Dirichlet boundary conditions. Here 
Bd(t) = Y^=i n d =i( n i + • • • + n d)~* is the generalized Epstein zeta function [96|. convergent for d < 2t. Of course, 
in a finite trap at temperatures of the order of or less than the energy of the first excited quasiparticle, T < ei, the 
scaling law ()294ll is no longer valid and the condensate fluctuations acquire a different temperature scaling, due to 
the temperature-independent quantum noise produced by the excited atoms, which are forced by the interaction to 
occupy the excited energy levels even at zero temperature, as it was discussed for the Eqs. (|279|l and (|244|l . 

4- Constraint mechanism of anomalous order-parameter fluctuations and susceptibilities 
versus instability in the systems with a broken continuous symmetry 

It is important to realize that the anomalously large order-parameter fluctuations and susceptibilities have a simple 
geometrical nature, related to the fact that the direction of the order parameter is only in a neutral, rather than in a 
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stable equilibrium, and does not violate an overall stability of the system with a broken continuous symmetry at any 
given temperature below phase transition, T < T c . On one hand, on the basis of a well-known relation between the 
longitudinal susceptibility and the variance Am of the order parameter fluctuations, 

1 r)M A 2 

Xaa = = J±M_ A 2 = ,t M _ m J 2 ), (296) 

A N dB a Nk B T' M u a a> h v ' 

it is obvious that if the fluctuations are anomalous, i.e., go as A M ~ N J with 7 > 1 instead of the standard macroscopic 
thermodynamic scaling A 2 ^ ~ N, then the longitudinal susceptibility diverges in the thermodynamic limit N — > 00 
as A\ [ /N ~ A^ 7 " 1 — > 00. On the other hand, the susceptibilities in stable systems should be finite, since otherwise 
any spontaneous perturbation will result in an infinite response and a transition to another phase. One could argue 
(as it was done in a recent series of papers [98] ) that the anomalous fluctuations cannot exist since they break the 
stability condition and make the system unstable. However, such an argument is not correct. 

First of all, the anomalously large transverse susceptibility in the infrared limit \± ~ k~ 2 ~ L 2 in Bogoliubov's 
l/fc 2 -theorem (|289ll originates from the obvious property of a system with a broken continuous symmetry that it is 
infinitesimally easy to change the direction of the order parameter and, of course, does not violate stability of the 
system. However, it implies anomalously large fluctuations in the direction of the order parameter. Therefore, an 
anomalously large variance of the fluctuations comes from the fluctuations (M — M) _L M, which are perpendicular to 
the mean value of the order parameter. This is the key issue. It means that a longitudinal external field B || M, that 
easily rotates these transverse fluctuations towards M when x±B ~ Am, will change the longitudinal order parameter 



by a large increment, x\\B ~ |M| — yM 2 — A 2 f A M /2M. It is obvious that this pure geometrical rotation of 
the order parameter has nothing to do with an instability of the system, but immediately reveals the anomalously 
large longitudinal susceptibility and relates its value to the anomalous transverse susceptibility and variance of the 
fluctuations, 

XII ~ A 2 M /2BM ~ AmxJZM. (297) 

This basically geometrical mechanism of an anomalous behavior of constrained systems constitutes an essence of 
the constraint mechanism of the infrared anomalies in fluctuations and susceptibilities of the order parameter for 
all systems with a broken continuous symmetry. The latter qualitative estimate yields the anomalous scaling of 
the longitudinal susceptibility discussed above, \\\ ~ L with increase of the system size L, Eq. (|291|l . since x± ~ 
L 2 , M ~ L 3 , and A M ~ L 2 . 

A different question is whether the Bose gas in a trap is unstable in the grand canonical ensemble when an 
actual exchange of atoms with a reservoir is allowed, and only the mean number of atoms in the trap is fixed, 
N = const. In this case, for example, the isothermal compressibility is determined by the variance of the number- 
of-atoms fluctuations kt = — V~ 1 (dV/dP)x — V{(N — N) 2 ) / N 2 ksT, and diverges in the thermodynamic limit if 
fluctuations are anomalous. In particular, the ideal Bose gas in the grand canonical ensemble does not have a well- 



defined condensate order parameter, since the variance is of the order of the mean value, y ((N — N) 2 ) ~ N, and it 
is unstable against a collapse 0, 0, |9?| . 

In summary, one could naively expect that the order parameter fluctuations below T c are just like that of a 
standard thermodynamic variable, because there is a finite restoring force for deviations from the equilibrium value. 
However, in all systems with a broken continuous symmetry the universal existence of infrared singularities in the 
variance and higher moments ensures anomalously large and non-Gaussian fluctuations of the order parameter. This 
effect is related to the long- wavelength phase fluctuations and the infrared singularity of the longitudinal susceptibility 
originating from the inevitable geometrical coupling between longitudinal and transverse order parameter fluctuations 
in constrained systems. 



VII. CONCLUSIONS 



It is interesting to note that the first results for the average and variance of occupation numbers in the ideal Bose 
gas in the canonical ensemble were obtained about fifty years ago by standard statistical methods [42ll75ll7r3 | (see also 
|74l l77| and the review IT^|). Only later, in the 60s, laser physics and its byproduct, the master equation approach, 
was developed (see, e.g.. |l6l K^). In this paper we have shown that the latter approach provides very simple and 
effective tools to calculate statistical properties of an ideal Bose gas in contact with a thermal reservoir. In particular, 
the results (|169fl and (|172|) reduce to the mentioned old results in the "condensed region" in the thermodynamic limit. 

However, the master equation approach gives even more. It yields simple analytical expressions for the distribution 
function of the number of condensed atoms (|162|l and for the canonical partition function i|163fl . In terms of cumulants, 
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or semi- invariants |72|, |95j, for the stochastic variables no or n — N — no, it was shown [2l| that the quasithermal 
approximation l|154l) . with the results (|169[) and l|172|l . gives correctly both the first and the second cumulants. The 
analysis of the higher-order cumulants is more complicated and includes, in principle, a comparison with more accurate 
calculations of the conditioned average number of non-condensed atoms (I129[l as well as higher-order corrections to 
the second order master equation l|124l> . It is clear that the master equation approach is capable of giving the 
correct answer for higher-order cumulants and, therefore, moments of the condensate fluctuations. Even without 
these complications, the approximate result l|162[l reproduces the higher moments, calculated numerically via the 
exact recursion relation (JSDJl, remarkably well for all temperatures T < T c and T ~ T c (see Fig. I12|) . 

As we demonstrated in Section IV, the simple formulas yielded by the master equation approach allow us to 
study mesosco pic effects in BECs for a relatively small number of atoms that is typical for recent experiments 
|24l |25L 12a. I27L 123 . |29|. Moreover, it is interesting in the study of the dynamics of BEC. This technique for studying 
statistics and dynamics of BEC shows surprisingly good results even within the sim plest ap proximations. Thus, the 
analogy with phase transitions and quantum fluctuations in lasers (see, e.g., [T^. I^TTol Flf ) clarifies some problems 
in BEC. The equilibrium properties of the number-of-condensed-atom statistics in the ideal Bose gas are relatively 
insensitive to the details of the model. The origin of dynamical and coherent properties of the evaporatively cooling 
gas with an interatomic interaction is conceptually different from that in the present "ideal gas + thermal reservoir" 
model. The present model is rather close to the dilute -He 4 gas in porous gel experiments |22| in which phonons in 
the gel play the role of the external thermal reservoir. Nevertheless, the non-condensed atoms always play a part of 
some internal reservoir, and the condensate master equation probably contains terms similar to those in Eq. I|13U|) 
for any cooling mechanism. 

For the ideal Bose gas in the canonical ensemble the statistics of the condensate fluctuations below T c in the 
thermodynamic limit is essentially the statistics of the sum of the non-condensed modes of a trap, n^, that fluctuate 
independently, fiQ — N ~ Sk^o" -14 - This is well understood, especially, due to the Maxwell's demon ensemble 
approximation elaborated in a series of papers 0, E2j 0, EE 0, E3, and is completed and justified to a certain 
extent in [2E by the explicit calculation of the moments (cumulants) of all orders, and by the reformulation of 
the canonical-ensemble problem in the properly reduced subspace of the original many-particle Fock space. The main 
result (|263|) of [20I l2ll| explicitly describes the non-Gaussian properties and the crossover between the ideal-gas and 
interaction-dominated regimes of the BEC fluctuations. 

The problem of dynamics and fluctuations of BEC for the interacting gas is much more involved. The master 
equation approach provides a powerful tool for the solution of this problem as well. Of course, to take into account 
higher-order effects of interaction between atoms we have to go beyond the second order master equation, i.e., to 
iterate Eq. 1123|) more times and to proceed with the higher-order master equation similarly to what we discussed 
above. It would be interesting to show that the master equation approach could take into account al l higher order 
effects in a way generalizing the well-known nonequilibrium Keldysh diagram technique |30l Il06t Il07j . As a result, 
the second order master equation analysis presented above can be justified rigorously, and higher-order effects in 
condensate fluctuations at equilibrium, as well as nonequilibrium stages of cooling of both ideal and interacting Bose 
gases can be calculated. 

The canonical-ensemble quasiparticle method, i.e., the reformulation of the problem in terms of the proper canonical- 
ensemble quasiparticles, gives even more. Namely, it opens a way to an effective solution of the canonical ensemble 
problems for the statistics and nonequilibrium dynamics of the BEC in the interacting gas as well. The first step in 
this direction is done in Section VI, where the effect of the Bogoliubov coupling between excited atoms due to a weak 
interaction on the statistics of the fluctuations of the number of ground-state atoms in the canonical ensemble was 
analytically calculated for the moments (cumulants) of all orders. In this case, the BEC statistics is essentially the 
statistics of the sum of the dressed quasiparticles that fluctuate independently. In particular, a suppression of the 
condensate fluctuations at the moderate temperatures and their enhancement at very low temperatures immediately 
follow from this picture. 

There is also the problem of the BEC statistics in the microcanonical ensemble, which is closely related to the 
canonical-ensemble problem. In particular, the equilibrium microcanonical statistics can be calculated from the 
canonical statistics by means of an inversion of a kind of Laplace transformation from the temperature to the energy 
as independent variable. Some results concerning the BEC statistics in the microcanonical ensemble for the ideal 
Bose gas were presented in 0, EH HE Eil EE EH ■ We can calculate all moments of the microcanonical fluctuations of 
the condensate from the canonical moments found in the present paper. Calculation of the microcanonical statistics 
starting from the grand canonical ensemble and applying a saddle-point method twice, first, to obtain the canonical 
statistics and, then, to get the microcanonical statistics |45| . meets certain difficulties since the standard saddle-point 
approximation is not always good and explicit to restore the canonical statistics from the grand canonical one with 
sufficient accuracy 57] . The variant of the saddle-point method discussed in Appendix IF1 is not subject to these 
restrictions. 

Another important problem is the study of mesoscopic effects due to a relatively small number of trapped atoms 
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(N ~ 10 3 — 10 6 ). The canonical-ensemble quasiparticle approach under the approximation (|210|l . i.e., 7i CE w W^Lg, 
takes into account only a finite-size effect of the discreteness of the sing le p article spectrum, but does not include all 
mesoscopic effects. Hence, other methods should be used (see, e.g., |l4l4l. l46ll4^ l6^l6^.l63.l63.l6^|). In particular, 
the master equation approach provides amazingly good results in the study of mesoscopic effects, as was demonstrated 
recently in [53, |53| • 

The canonical-ensemble quasiparticle approach also makes it clear how to extend the Bogoliubov and more advanced 
diagram methods for the solution of the canonical-ensemble BEC problems and ensure conservation of the number 
of particles. The latter fact cancels the main arguments of Refs. |97J, |98j a gain st the Giorgini-Pitaevskii-Stringari 
result 78] and shows that our result l|263|) and, in particular, the result of |78j ] for the variance of the number of 
ground-state atoms in the dilute weakly interacting Bose gas, correctly take into account one of the main effects of 
the interaction, namely, dressing of the excited atoms by the macroscopic condensate via the Bogoliubov coupling. If 
one ignores this and other correlation effects, as it was done in [13,113, the result cannot be correct. This explains 
a sharp disagreement of the ground-state occupation variance suggested in [92J with the predictions of [78| and our 
results as well. Note also that the statement from that "the phonon spectrum plays a crucial role in the approach 
of [Z^" should not be taken literally since the relative weights of bare modes in the eigenmodes (quasiparticles) is, 
at least, no less important than eigenenergies themselves. In other words, our derivation of Eq. (|263[1 shows that 
squeezing of the excited states due to Bogoliubov coupling in the field of the macroscopic condensate is crucial for 
the correct calculation of the BEC fluctuations. Besides, the general conclusion that very long wavelength excitations 
have an acoustic, "gapless" spectrum (in the thermodynamic limit) is a cornerstone fact of the many-body theory of 
superfluidity and BEC 94] . Contrary to a pessimistic picture of a mess in the study of the condensate fluctuations 
in the interacting gas presented in |97|. we are convinced that the problem can be clearly formulated and solved 
by a comparative analysis of the contributions of the main effects of the interaction in the tradition of many-body 
theory. In particular, the result (|263l) corresponds to a well established first-order Popov approximation in the diagram 
technique for the condensed phase |36j . 

We emphasize here an important result of an analytical calculation of all higher cumulants (moments) |2l| • In most 
cases (except, e.g., for the ideal gas in the harmonic trap and similar high dimensional traps where d > 2a), both for 
the ideal Bose gas and for the interacting Bose gas, the third and higher cumulants of the number-of-condensed-atom 
fluctuations normalized to the corresponding power of the variance do not tend to the Gaussian zero value in the 
thermodynamic limit, e.g., ((no — (no)) 3 ) /((no — ( n o)) 2 ) 3 ^ 2 does not vanish in the thermodynamic limit. 

Thus, fluctuations in BEC are not Gaussian, contrary to what is usually assumed following the Einstein theory of 
fluctuations in the macroscopic thermodynamics. Moreover, BEC fluctuations are, in fact, anomalously large, i.e., 
they are not normal at all. Both these remarkable features originate from the universal infrared anomalies in the 
order-parameter fluctuations and susceptibilities in constrained systems with a broken continuous symmetry. The 
infrared anomalies come from a long range order in the phases below the critical temperature of a second order phase 
transition and have a clear geometrical nature, related to the fact that the direction of the order parameter is only in a 
neutral, rather than in a stable equilibrium. Hence, the transverse susceptibility and fluctuations are anomalously large 
and, through an inevitable geometrical coupling between longitudinal and transverse order parameter fluctuations in 
constrained systems, produce the anomalous order parameter fluctuations. In other words, the long wavelength phase 
fluctuations of the Goldstone modes, in accordance with the Bogoliubov l/fc 2 -theorem for the transverse susceptibility, 
generate anomalous longitudinal fluctuations in the order parameter of the systems below the critical temperature of 
the second order phase transition. Obviously, this constraint mechanism of the infrared anomalies in fluctuations and 
susceptibilities of the order parameter is universal for all systems with a broken continuous symmetry, including BEC 
in ideal or weakly interacting gases as well as superfluids, ferromagnets and other systems with strong interaction. 
It would be interesting to extend the analysis of the order parameter fluctuations presented in this review from the 
BEC in gases to other systems. 

The next step should be an inclusion of the effects of a finite renormalization of the energy spectrum as well as the 
interaction of the canonical-ensemble quasiparticles at finite temperatures on the statistics and dynamics of BEC. It 
can be done on the level of the second-order Beliaev-Popov approximation, which is considered to be enough for the 
detailed account of most many-body effects (for a review, see |36j). A particularly interesti ng p roblem is the analysis of 
phase fluctuations of the condensate in the trap, or of the matter beam in the atom las er 1891, because the interaction 
is crucial for the existence of the coherence in the condensate [2I M 153, M MM fTosTlinfll [Tirl| . As far as the 
equilibrium or quasi-equilibrium properties are concerned, the problem can be solved effectively by applying either the 
traditional methods of statistical physics to the canonical-ensemble quasiparticles, or the master equation approach, 
that works surprisingly well even without any explicit reduction of the many-particle Hilbert space [5^ . l53| . For the 
dynamical, nonequilibrium properties, the analysis can be bas ed on an appropriate modification of the well-known 
nonequilibrium Keldysh diagram technique flM Hoi HT1 HT^ which incorporates both the standard statistical and 
master equation methods. 

Work in the directions mentioned above is in progress and will be presented elsewhere. Clearly, the condensate and 
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non-condensate fluctuations are crucially important for the process of the second order phase transition, and for the 
overall physics of the Bose-Einstein-condensed interacting gas as a many-particle system. 
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APPENDIX A: BOSE'S AND EINSTEIN'S WAY OF COUNTING MICROSTATES 

When discussing Einstein's 1925 paper [jj, we referred to the identity (|49|) . 

Z\ {N + Z-l 



p \ ... p N \ v N 

{P0,P1,—,PN} 



(Al) 



which expresses the number of ways to distribute N Bose particles over Z quantum cells in two different manners: 
On the left-hand side, which corresponds to Bose's way of counting microstates, numbers p r specify how many cells 
contain r quanta; of course, with only N quanta being available, one has p r = for r > N. As symbolically indicated 
by the prime on the summation sign, the sum thus extends only over those sets {po,Pi, ■ ■ • ,Pn} with comply with 
the conditions 



Y,Pr = Z, (A2) 



r=0 

stating that there are Z cells to accommodate the quanta, and 

N 



r Pr = N , (A3) 



E 

r=0 

stating that the number of quanta be N. The right-hand side of Eq. I|A1|) gives the total number of microstates, 
taking into account all possible sets of occupation numbers in the single expression already used by Einstein in his 
final derivation Q of the Bose-Einstein distribution which can still be found in today's textbooks. 

The validity of Eq. (|A1I) is clear for combinatorial reasons. Nonetheless, since this identity i|Al(l constitutes one of 
the less known relations in the theory of Bose-Einstein statistics, we give its explicit proof in this appendix. 

As for most mathematical proofs, one needs tools, an idea, and a conjurer's trick. In the present case, the tools are 
two generalizations of the binomial theorem 

{a + b) n = Y. (f)a k b n - k , (A4) 



.k 

k=Q v 



where 



(A5) 



y kJ k\(n-k)\ 

denotes the familiar binomial coefficents. The first such generalization is the multinomial theorem 

(a 1+ a 2 + ... + a N ) n = V — — ^ -afaf...<, (A6) 

E Vr=TL 

which is easily understood: When multiplying out the left-hand side, every product obtained contains one factor 
dj from each bracket (a\ + a 2 + . . . + ajv). Hence, in every product the exponents add up to the 

number of brackets, which is n. Therefore, for such a product there are n\ permutations of the individual factors a\. 
However, if identical factors di are permuted among themselves, for which there are pi\ possibilities, one obtains the 
same value. Hence, the coefficient of each product on the right-hand side in Eq. l|A6jl corresponds to the number of 
possible arrangements of its factors, divided by the number of equivalent arrangements. Note that the reasoning here 
is essentially the same as for the justification of Bose's expression <|12[) . which is, of course, not accidental. 
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The second generalization of the binomial theorem (|A4|) required for the proof of the identity ijAlfl emerges when 
we replace the exponent n by a non-natural number 7: One then has 

(a + &) 7 = £(I) a * 67 " fc ' (A7) 

with the definition 

f i\ = 7(7-l)(7-2)-...-( 7 -fc + l) _ (Ag) 
fey fe! 

If 7 is not a natural number, this series (|A7() converges for any complex numbers a, b, provided \a/b\ < 1. This 
generalized binomial theorem ()A7|) . which is treated in introductory analysis courses, is useful, e.g., for writing down 
the Taylor expansion of (1 + x) 1 . 

Given these tools, the idea for proving the identity IA1I now consists in considering the expression 

Pnz{x) = (1 + x + x 2 + x 3 + . . . + x N ) z , (A9) 

where x is some variable which obeys |x| < 1, but need not be specified further. According to the multinomial 
theorem (|A6|I . one has 

Pnz(x)= V -j— ^ :(x°Y»(x l Y* ...(x N )™ . (A10) 

Directing the attention then to a systematic ordering of this series with respect to powers of x, the coefficient of x N 
equals the sum of all coefficients encountered here which accompany terms with ^ r rp r = N, which is precisely the 
left-hand side of the desired equation (|A1|) . keeping in mind the restrictions i|A2(l and l|A3(l . 

Now comes the conjurer's trick: The coefficient of x N in this expression (|A10(I . equaling the left-hand side of 
Eq. IjAljl . also equals the coefficent of x N in the expression 

Poo(x) = (1 + x + x 2 + x 3 + . . .x N + x N+1 + . . .) z , (All) 

since this differs from Pnz{%) only by powers of x higher than N . But this involves a geometric series, which is 
immediately summed: 

/ ~, \ z 



Poo(x) = E : 

\r=0 / 

= (1 -x)- z . (A12) 
According to the generalized binomial theorem (| AT|> . one has 

a -*r*=f; ("5) (-*)*. ( ai3 ) 

N 



so that, in view of the definition (|A8|) . the coefficient of x equals 



(-l)^(-^ = { _ l)N (-Z)(-Z-l)-...-{-Z-N + l) 



N J Nl 

{Z + N-l)- ...-(Z + 1)Z 
N\ 

(Z + N- 1)! 
TV! (Z-l)! 

This is the right-hand side of the identity I |A1| I. which completes the proof. 
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APPENDIX B: ANALYTICAL EXPRESSION FOR THE MEAN NUMBER OF CONDENSED ATOMS 

We obtain an analytical expression for no from Eq. (|86() . One can reduce the triple sum into a single sum if we take 
into account the degeneracy g(E) of the level with energy E = Ml(l + m + n), which is equal to the number of ways 
to fill the level. This number can be calculated from the Einstein's complexion equation, i.e. 

{(Z + rn + n) + (3-1)}! l(E \(E \ 
9{E) = (Z + m + n)!(3-l)l = 2^ + jp + T ( } 



Now, we have reduced three variables /, m, n to only one. By letting E = shfl where s is integer, one can write 

f^ Q f± + l) e »/J«i _ i (1 + no) 

where for no 3> 1 



i( S + 2)( S + l) _ n 

n + _ S, (B2) 



S— 1 



The root of the quadratic Eq. (|B2I) yields 



no = -^[(1 + S - N) - V(l + S - N) 2 + AN}. (B4) 



In order to find an analytical expression for S, we write 

00 D 2 



lv^ S + OS 1 

s=0 s=l 



Converting the summation into integration by replacing a; = s/3Ml yields 

I f°° dx 1 f°° x 2 dx 3 



1 2a 3 ,/ e x ~ 1 2a 



- l-^-D + ^+B 2 , (B6) 



a a 3 V 2a 

where a = j3Ml = Kl/k B T = (((3)/N) 1 / 3 T C /T. 

Thus, Eq. I|B3|I gives the following analytical expression for S: 

(1 + S-N) 



-Nil- 



T\ 3 1 7r 2 /iV\ 2 ^ 3 /T x2 



Tc) j 4 vc(3)y \T C 

V /d ^ln ( e («3)/iV)^T c /T_^ +2 



U(3), 

Fig. 1161 compares different approximations for calculating h within the grand canonical ensemble. 



(B7) 



APPENDIX C: FORMULAS FOR THE CENTRAL MOMENTS OF CONDENSATE FLUCTUATIONS 

By using no — N — Y n k and (no) = N — Y ( n k), we have 
fe^o fc/0 

{(no (no)) s ) = (-l)'(EK - (n fe ))H, (CI) 

k^O 
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FIG. 16: Grand canonical result for (no) as a function of temperature for N = 200, computed from analytical expressions, Eqs. 

and (1 ine with circles); semi-analytical expressions given by Eqs. lB4t and 1B3II (solid line); and exact numerical 

solution of Eq. (1861 (dots) . Expanded views show: b) exact agreement of the semi-analytical approach at low temperature and 
c) a small deviation near T c . Also, the analytical and the semi-analytical results agree quite well. 



which shows that the fluctuations of the condensate particles are proportional to the fluctuations of the non-condensate 
particles. 

As an example we show how to evaluate the fluctuation for the second order moment, or variance, 



((n - (n )) 2 } = (( n jn k ) ~ (n 3 )(n k )). 



(C2) 



The numbers of particles in different levels are statistically independent, since (rijUk^j) = Tr{dMjd k dkp} = 



Tr{alajPj}Tr{al<ikpk} = (nj)(n k ). Thus, we find 



((«o - 


(no)) 2 ) 




(no)) 3 ) 




(no)) 4 ) 



k^O 



k^O 



k^O 



In the grand canonical approach, (n k s ) can be evaluated using 



\ n k ) = tt > n h e 



(C3) 
(C4) 
(C5) 

(C6) 



where Z k = (1 — e /3 ^ £k ^) 1 . An alternative way is to use the formula (n k s ) = %(®y \u=o derived in Section V.B. In 
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particular, one can show that in this approach 



(n 2 k ) = 2(n k ) 2 + (n k ), 

(n|) = 6(n fc ) 3 + 6(n fc ) 2 + (n fc ), 

(4) = 24(n fe ) 4 + 36(n fe ) 3 + 14(n fc > 2 + (re- 



using Eqs. lfn7) - lfn9|l we obtain 

((n - (n )) 2 



fc^O 

((n -(n )) 3 ) = -]T[2(n fc ) 3 + 3(n fc ) 2 + (n fc )] 
kjto 

((no - («o)) 4 ) = E[ 9 ^ fe ) 4 + 18 < nfc ) 3 + 10 < rifc ) 2 + ^ fe )] 



(C7) 
(C8) 
(C9) 



(CIO) 



(Gil) 



APPENDIX D: ANALYTICAL EXPRESSION FOR THE VARIANCE OF CONDENSATE 

FLUCTUATIONS 

In a spherically symmetric harmonic trap with trap frequency fl one can convert the triple sums in Eq. I|98|l into a 

oo 

single sum by using Eq. (|B1J| . and then do integration upon replacing Y] ... — > J Q dx... and (3E —* x, giving 



IE=0 



£#0 



2 V nn 



Ml 



e X p(/3^)(l + J-)-l 



ex P (/3I?) 1 + i -1 



1 



dx 



( x 



■ix 



\pnnj phfi 



exp(x) i + J- -1 



+ 



exp 



0) (l + - 



- 1 







3a; 


— / dx 


[(;)"< 


- — + 2 


2« A 




a 



exp 



exp^) (l + 



(Dl) 



where a = (3MI = hQ/k B T = (T C /T)(C(3)/7V) 1/3 and the density of states can be written as p{E) 
. Then we integrate by parts using the identity 



l 

2M2 



'E_y + M + 9 

v Tin / ' an ' 



exp 



(x) f 1 + J- 



exp 



- 1 



c9x 



exp(x) (l + JL -1 



arriving at 



An 



a 


e a (l ^ 




- 1 


2a 2 




"0 J 







dx ■ 



+3) 



exp 



- 1 



The integral in Eq. (|D2|> can be calculated analytically, using 

xdx it 2 1 



a [exp(x)A-l] 6 2 



In 2 A + In ( Ae a - 1) hi A + di log(^e a ) 



(D2) 



72 



where 



As a result, we get 



Ang = \ 
a 6 



f°° dx 

I leM*)A-l} =Hl+a) - MAea - 1)+a > 

ln(i) 



di log(z) 



di log 



e a 1 



In" 1 



1-t 



n 



-dt. 



In 



e a 1 



"o 



In 1 



3 . / n + l 
-am — r 

2 \e a (no + l)-n 



a 


e a (l 


+ &) 


- 1 











Taking into account that a = (T C /T)(C(3)/7V) 1 / 3 , we finally obtain 



Arc* = [ L 



N 



TcJ C(3) 



— + di log 
6 



exp[(T c /T)(C(3)/7V) 1 /3] (l + L 
V no 



2 V «o 



(D3) 



In 



exp[ 



( 1+ n J 


- 1 











?l + 1 



2 VT c y VC(3)y Lexp[(T c /T)(C(3)/iV)V3](n + l)-no 



T f N 



T c VC(3) 



1/3 



exp[(T c /T)(C(3)/iV)V3] (l + ^) - 1 



(D4) 



APPENDIX E: SINGLE MODE COUPLED TO A RESERVOIR OF OSCILLATORS 



The derivation of the damping Liouvillean proceeds from the Liouville-von Neumann equation 

Ft m = Jh [fsr ' mh (E1) 

where V sr is the Hamiltonian in the interaction picture for the system (s) coupled to a reservoir (r). 

We can derive a closed form of the dynamical equation for the reduced density operator for the system, 
p s (t) = Tr r {/5(i)} by tracing out the reservoir. This is accomplished by first integrating Eq. (|E1|) for p(t) = 
TE Jo p(t')]dt' and then substituting it back into the right hand side of Eq. IIElfl . giving 

Zps(t) = ^Tr r [V sr ,P(0)] + T-^Tr r f [V sr (t), [V sr (t'), p(t')]dt'] . (E2) 
at in (thy j Q 

We may repeat this indefinitely, but owing to the weaknesses of the system-reservoir interaction, it is possible to 
ignore terms higher than 2nd order in V sr . Furthermore, we assume the system and reservoir are approximately 
uncorrelated in the past and the reservoir is so large that it remains practically in thermal equilibrium p l r h , so pit') ~ 
Ps {t')® pf and p{0) = Ps (0)® pf. 

For a single mode field (/) coupled to a reservoir of oscillators, one has 

Vfr = K^9±{Ka ] e i{v ~ Vk)t + blae-*^-"^) . (E3) 

k 
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Since Tr r {a^/5*' 1 } = and Tr r {a/5' fc } = 0, the first term vanishes and by using Eq. I|E3|) we have 16 terms. But secular 
approximation reduces the number of terms by half. We now perform the Markov approximation, pf(t') ~ Pf(t), 
stating that the dynamics of the system is independent of the states in the past. 

The thermal average of the radiation operators is 

k J ° k 

where Tr r {pfb1J)v,} = n(z^) = {e 0hVk - l)" 1 . Thus, we have 

§- t Pf(t) = -\c[rfap f {i) - 2ap f (t)tf + Pf(t)a!a] - ±D[atfp f {t) - 2atp f (t)a + p f (t)atf] , (E4) 
where T> = Qn and C = G(n + 1). 



APPENDIX F: THE SADDLE-POINT METHOD FOR CONDENSED BOSE GASES 



The saddle-point method is one of the most essential tools in statistical physics. Yet, the conventional form of this 
approximation fails in the case of condensed ideal Bose gases [Hi3. The point is that in the condensate regime 
the saddle-point of the grand canonical partition function approaches the ground-state singularity at z = exp(/3eo), 
which is a hallmark of BEC. However, the customary Gaussian approximation requires that intervals around the 
saddle-point stay clear of singularities. Following the original suggestion by Dingle [59f. Holthaus and Kalinowski |(3(j 
worked out a natural solution to this problem: One should exempt the ground-state factor of the grand canonical 
partition function from the Gaussian expansion and treat that factor exactly, but proceed as usual otherwise. The 
success of this refined saddle-point method hinges on the fact that the emerging integrals with singular integrands can 
be done exactly; they lead directly to parabolic cylinder functions. Here we discuss the refined saddle-point method 
in some detail. 

We start from the grand canonical partition function 

oo 1 

^ ,|= n .-H-w (F1) 

where e v are single-particle energies, (3 — l/ksT and z — exp(/3/i). The grand canonical partition function 5(/3, z) 
generates the canonical partition functions (/?) by means of the expansion 

oo 

SO?,*) = 5>*%r(/3). (F2) 

Then we treat z as a complex variable and using Cauchy's theorem represent Zn(0) by a contour integral, 

ZN{f3) ^ in f dZ ^^ ^ ^-fdzeM-F(z)), (F3) 
where the path of integration encircles the origin counter-clockwise, and 

oo 

F(z) = (iV + l)ln^-ln5(/3,z) = (N + 1) lnz + ^ M 1 ~ zexp(-0e v )]. (F4) 

The saddle-point zo is determined by the requirement that F(z) becomes stationary, 

dF(z) 



dz 



= 0, (F5) 

z=z 



giving 



2V + 1 = Y] j--— t -. (F6) 



v=0 
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FIG. 17: Canonical occupation number of the ground state as a function of temperature for N = 200 atoms in a harmonic 
isotropic trap. Dashed and solid curves are obtained by the conventional and the refined saddle-point method, respectively. 
Dots are "exact" numerical answers obtained for the canonical ensemble. 



This is just the grand canonical relation between particle number N and fugacity zq (apart from the appearance of 
one extra particle on the left-hand side). Appendix C discusses an approximate analytic solution of such equations. 
In the conventional saddle-point method, the whole function F(z) is taken in the saddle-point approximation 

F(z)*F(z ) + If"(z )(z-z )2. (F7) 
Doing the remaining Gaussian integral yields 

Z N {(3) « e " p( ~ F(Zo)) . (F8) 
V-2tt^"(zo) 

The canonical occupation number of the ground state, and its mean-square fluctuations, are obtained by differentiating 
the canonical partition function: 

_ _ d\nZ N (p) _ 1 1 r 1 05(/?,z) 

no ~ d(-/3s ) ~ Z N {f3)2mf aZ z^d{-l3E Q y [ " J) 



The saddle-point approximation is then applied to the integrands of Eqs. (|F9|) and (|F10|) . Figs. 1171 and 1181 (dashed 
curves) show results for no and Ano obtained by the conventional saddle-point method for a Bose gas with N = 200 
atoms in a harmonic isotropic trap. In the condensate regime there is a substantial deviation of the saddle-point 
curves from the "exact" numerical answer obtained by solution of the recursion equations for the canonical statistics 
(dots). 

The reason for this inaccuracy is that in the condensate region the saddle-point zq lies close to the singular point 
z = exp(/3eo) of the function F(z). As a result, the approximation i|F7|l becomes invalid in the condensate region. 
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FIG. 18: Variance in the condensate particle number as a function of temperature for N = 200 atoms in a harmonic isotropic 
trap. Dashed and solid curves are obtained by the conventional and the refined saddle-point method, respectively. Dots are 
"exact" numerical answers in the canonical ensemble. 



To improve the method, Dingle [S9| proposed to treat the potentially dangerous term in (|F4|) as it is, and represent 
Z N (P) as 

^) = ^Uz-^tIML (F11) 

2m J 1 — zexp(— (jeo) 

where 

oo 

Fi(z) = (iV+l)lnz + ^ln[l-zexp(-/3e t/ )] (F12) 

!/=l 

has no singularity at z — exp(/3eo)- The singular point to be watched now is the one at z — exp(/3ei). Since 
z < exp(/3e ), the saddle-point remains separated from that singularity by at least the TV-independent gap exp(/3ei) — 
exp(/3e ) — huj/ksT. This guarantees that the amputated function Fi(z) remains singularity-free in the required 
interval around zq for sufficiently large N. Then the Gaussian approximation to exp(— F\(z)) is safe. The subsequently 
emerging saddle-point integral for the canonical partition function can be done exactly, yielding [6fj| | 

Z N (J3) « i exp[/3e - F x (z ) - 1 + r; 2 /2]erfc (^j , (F13) 



where erfc(z) = 2/a/tt J z °° exp(—t 2 )dt is the complementary error function, r/ = (exp(/3eo) — z a) \/ ~ F" ( z o) , an d 
T}i=V- 1 /V- 

Calculation of occupation numbers and their fluctuations deals with integrals from derivatives of S(/3, z) with respect 
to —(3sq. In such expressions the factors singular at z = exp(/3eo) should be taken exactly. This leads to the integrals 
of the following form: 

h j ^ ^feipt-w' 01 » jgl-KMr-™ «p(*o -iW-* n'ft - (F14) 
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FIG. 19: The third central moment {(no — no) 3 } for N — 200 atoms in a harmonic isotropic trap obtained using the refined 
saddle-point method (solid curve). Exact numerical results for the canonical ensemble are shown by dots. 



where i] = (exp(/3eo) ~ z *)\/~ f'{{ z *)i Vi = V ~ a Mi an d is a saddle-point of the function 

f(z) = h(z) + (<r- l)0e o +<rln(l -zexp(-0e o )); 
D_ CT (z) is a parabolic cylinder function, which can be expressed in terms of hypergeometric functions as 

'iFiC-sA 1/2, z 2 /2) V2z •! Fi((l - s)/2, 3/2, z 2 /2)~ 



_ 9 s/2 -z 2 /4 



na 



0/2] 



r[-.s/2] 



(F15) 



(F16) 



Figs. 1171 and 1181 (solid curves) show no(T) and Ano(T) obtained by the refined saddle-point method. These results 
are in remarkable agreement with the exact dots. Fig. 1191 demonstrates that this refined method also provides good 
accuracy for the third central moment of the number-of-condensed-atoms fluctuations. 
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